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ABSTRACT 

A three-dimensional MHD model for the propagation and dissipation of Alfven waves in a coronal loop is 
developed. The model includes the lower atmospheres at the two ends of the loop. The waves originate on small 
spatial scales (less than 100 km) inside the kilogauss flux elements in the photosphere. The model describes the 
nonlinear interactions between Alfven waves using the reduced MHD approximation. The increase of Alfven 
speed with height in the chromosphere and transition region (TR) causes strong wave reflection, which leads to 
counter-propagating waves and turbulence in the photospheric and chromospheric parts of the flux tube. Part 
of the wave energy is transmitted through the TR and produces turbulence in the corona. We find that the hot 
coronal loops typically found in active regions can be explained in terms of Alfven wave turbulence, provided 
the small-scale footpoint motions have velocities of 1-2 km/s and time scales of 60-200 s. The heating rate per 
unit volume in the chromosphere is 2 to 3 orders of magnitude larger than that in the corona. We construct a 
series of models with different values of the model parameters, and find that the coronal heating rate increases 
with coronal field strength and decreases with loop length. We conclude that coronal loops and the underlying 
chromosphere may both be heated by Alfvenic turbulence. 

Subject headings: Sun: atmospheric motions — Sun: chromosphere — Sun: corona — Sun: magnetic fields 
— turbulence — magnetohydrodynamics (MHD) 



1. INTRODUCTION 

It has long been assumed that the solar corona is heated by 
dissipation of magnetic disturba nces that prop agate up from 
the Sun's convection zone (e.g.. lAlfvenlll94'7h . Convective 
flows interacting with magnetic flux elements in the photo- 
sphere can produce magneto-hydrodynamic (MHD) waves 
that propagate up along the flux tubes and dissipate their en- 
ergy in the corona. Also, in closed magnetic structures the 
random motions of photospheric footpoints can lead to twist- 
ing and braiding of coronal field lines, an d to the formation of 
thin current sheets in the corona (also see lParkerilT972l [T9831 
119941: IPriest et al.ll2002h . Magnetic reconnection within such 
current sheet s may cause impulsive heating events, called 
"nanoflares" (P arker! [l988l) . The observed X-ray emission 
from solar active regions may be due to the cumulative ef- 
fects of many such coronal heating events. However, the de- 
tailed physical processes by which the corona is heated are 
not yet fully understood (for review s of coronal heating, see 
lAschwande n 2005':' Klimchu3l2006l) . The heating of solar ac- 
tive regions has in principle two contributions: (1) energy may 
be injected into the corona as a result of small-scale, random 
footpoint motions, or (2) the dissipated energy may originate 
from a large- scale nonpotential magnetic field suc h as a coro- 
nal flux rope (van Balleg ooijen & Cranm er 2008). In the sec- 
ond case the magnetic free energy is already stored in the 
corona, and does not need to be transported into the corona 
from the lower atmosphere. In this paper we will focus on the 
first case, i.e., we assume that the large-scale magnetic field 
of the active region is close to a potential field, and that most 
of the energy for coronal heating is provided by small-scale 
footpoint motions. 

Detailed MHD models of magnetic braiding have be en de- 
veloped by many authors (e.g., van Ballegooijen 1985. 119861 
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Il994t IHendrix & van Hoven 
Galsgaard & Nordlund 1996; Ng & Bhattachariee 1998 
Craig & Snevd 2 0051; iGudiksen & Nordlundi i2005 
Rappazzo et al. 2007n200l7~ This modeling has gener- 
ally confirmed Parker's ideas concerning the development 
of thin current layers in magnetic fields subject to ran- 
dom footpoint motion s . Ho wever, except for the work by 
iGudiksen & Nordlundi (l2005l) . the above-mentioned models 
do not attempt to describe the real structure and dynamics 
of the photospheric magnetic field. Most models use a 
highly simplified geometry in which the curvature of the 
coronal loops is neglected and the background magnetic 
field is assumed to be uniform, Bq = Bqz, where z is the unit 
vector in a Cartesian reference frame (Parker 1972). The 
photosphere at the two ends of the loop is represented by two 
parallel plates z = and z = L, and the magnetic field lines are 
assumed to be "line-tied" at these boundaries. The imposed 
horizontal flows (vi, v,.) at the boundaries are usually taken to 
be incompressible. This is very different from the converging 
and diverging motions observed at the solar surface. Obser- 
vations show that the photospheric magnetic field outside 
sunspots is highly intermittent and is concentrated into small 
flux elements or "flux tubes" with kilogauss field strengths 
and w idths of a few 100 km (Stenflo 1989; iBerger & Tid3 
1200 lb . These flux tubes are located in intergranular lanes and 
are continually jostled about by convective flows below the 
photosphere. These features of the boundary motions are not 
yet accurately represented in the braiding models. 

Another important aspect of the lower atmosphere is that 
it takes a significant amount of time for the effects of the 
footpoint motions to be transmitted into the corona. The 
flux tubes interact with convective flows at the base of the 
photosphere, and it takes 60 to 80 seconds for an Alfven 
wave to travel from that level to the base of the corona (a 
height difference of about 2 Mm). Furthermore, Alfven waves 
are s ubject to strong wave reflection (Ferraro & Pl umptonI 
Il958l) . as the Alfven speed increases from about 15 km/s in 
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the photospheric flux tubes to more than 1000 km/s in the 
corona. On first impact with the chromosphere-corona tran- 
sition region (TR), only a small fraction of the wave energy 
is actually transmitted into the co rona (e.g., Hollweg 1981|; 
ICranmer & van B allegooijen 2005). so it takes many bounces 
for the waves to be transmitted. Therefore, the relationship 
between the photospheric footpoint motions and the horizon- 
tal motions in the low corona is very complex and is subject 
to significant time delays. These delays are not included in 
existing models of field-line braiding. 

In previous models of coronal magnetic braiding it was 
assumed that the photospheric footpoint motions relevant to 
braiding occur on a horizontal length scale comparable to 
that of the solar granulation or larger (^^ > 1 Mm). Granula- 
tion flow patterns evolve on a time scale of a few minutes, 
and the kilogauss flux elements (located in the intergranu- 
lar lanes) are forced to move with these evolving flows. It 
was assumed that the main effect of the granulation is to pro- 
duce random displacements of the flux elements as a whole. 
However, there may also exist transverse motions inside the 
flux elements on a scale less than the element width (about 
100 km). This is plausible because the flux tubes are sur- 
rounded by convectiv e downflows that a r e highly turbulent 
fcattaneo et al. 200 31; IVogleret all 120051: IStein & Nordlundl 
|2()06; Bus hbv et al.ll2008h . In the region just below the pho- 
tosphere a flux element is subject to transverse motions that 
not only displace the flux element as a whole, but may also 
distort its shape and cause random intermixing of the field 
lines inside it. The process is illustrated in Figure [T] which 
shows a kilogauss flux tube being distorted by convective 
flows that push the field lines from left to right in the fig- 
ure. At present, little is known from observations about the 
magnitude of such transverse motions or the time scales in- 
volved. However, such small-scale transverse motions could 
have important effects on the upper atmosphere by produc- 
ing Alfven waves that propagate upward along the field lines 
and dissipate their energy in the chromosphere and corona. 
Alfven waves have i ndeed been ob served both in the chro- 
mosphere ( De Ponti eu et al.ll2007bl) and in the corona (e.g., 
iTomczyke t al.n2007i) . The purpose of the present paper is to 
investigate the possible role of small-scale random motions 
inside photospheric flux elements in heating the solar atmo- 
sphere. 

The chromosphere is a conduit for the transport of mass 
and energy into the corona. Actually, only a small frac- 
tion of the non-thermal energy injected into the solar atmo- 
sphere is transmitted to the corona; most of the energy is 
dissipated in the lower atmosphere. Therefore, to under- 
stand the heating of the Sun's upper atmosphere it is im- 
portant to study the structure, dynamics and heating of the 
chromosphere. High-resolution observations of the chromo- 
sphere have shown that it has a complex thermodynamic 
structure that is strongly in fluenced by the prese nce of mag- 
netic fields (see reviews by lJudgdIIOOet iRuttenlllOOTi) . The 
chromosphere is highly dynamic, and is filled with jet-like 
features such as mottles and dynam ic fibrils on the so lar disk 
(iRouppe v an der Voo rt et al. 2007; De Pontieu et al. 2007a') 
and spicules at the limb ( De Pontieu et al. 2007b) . Realis- 
tic three-dimensional MH D models of spicule-like stru ctures 
have been developed (e.g.. lMartinez-Svkoraet al.ll2009 ). In- 
ternetwork regions on the quiet Sun are affected by p-mode 
waves that travel upward from the photospher e and pro- 
duce s hocks that cause intermitte nt heating (Carlsson & SteinI 
Il997t lUlmschneider et al1l2005h . These shocks produce dis- 
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Fig. 1 . — Interaction of a magnetic flux element with convective flows in 
an intergranular lane. The red object indicates the magnetic element contain- 
ing a nearly vertical magnetic field, as indicated by the black an'ows. The 
blue arrows indicate the convective flow, which push on the flux tube from 
one side. Due to the stiffness of the magnetic field, horizontal momentum 
is transported upward, which results in a distortion of the shape of the flux 
tube and generates transverse motions inside it (green an'ows). We suggest 
these transverse motions create Alfven waves that propagate into the upper 
atmosphere. 

tinct asymmetries in the profiles of the Ca II H & K lines. 
However, the magnetic network and plage regions appear 
to be heated in a different way (Hasan & van Ballegooije^ 
l2008l) . First, these regions are continually bright in the cores 
of the Ca II H & K and Ca II 8542 A fines, and the width of the 
Ha li ne is enhanced co mpared to the non-magnetic surround- 
ings (Cauz zi et al.ll2009i) . indicating that the magnetic chro- 
mosphere is heated by a sustained heating process. Second, 
the wavelength profiles of the Ca II H line from network/plage 
elem ents are symmetric with respect to the rest wavelength 
(e.g.. lLites et al.ll993l:[Sheminova et"ani2 005). indicating that 
the heating is not due to acoustic shock waves. In this paper 
we investigate whether the heating of the magnetized chro- 
mosphere ma};_bedueJodissi£atron of Alfven waves as sug- 
gested bv lDe Pontieu et alj (l2007bl) . 

The paper is organized as follows. In Section 2 the obser- 
vational constraints on braiding models are discussed. It is 
shown that, if the corona is heated by dissipation of braided 
fields, the braiding must occur on small transverse length 
scales (less than a few arcseconds). This motivates us to de- 
velop a 3D MHD model for the dynamics of Alfven waves 
inside a magnetic flux tube extending from the photosphere 
through the chromosphere into the corona. The MHD model 
is described in Section 3, and modeling results are presented 
in Section 4. The results are further discussed in Section 5. 

2. OBSERVATIONAL CONSTRAINTS ON MAGNETIC BRAIDING 
MODELS 

2. 1 . Searching for Evidence of Braided Fields 

Active regions contain loop-like structures that are aligned 
with the direction of the coronal magnetic field. I n this sec- 
tion we use data from the X-Ray Telescope (XRT, Golu bet aP 
120071; iKano et aL 2008) on Hinode cKosugi et al, 2007h to 
search for braided magnetic fields in the corona. Figure 2i 
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Fig. 2.— (a) XRT image of active region AR 1 1041 on 2010 January 25, 
taken with Ti-poly filter. The image is an average over the period 11:02 UT 
to 1 1:52 UT, and shows the logaiithm of the intensity displayed as a negative 
(black and white differ a factor 100 in intensity). The .v and y coordinates are 
relative to solar disk center, (b) Spatially filtered version of the same image 
(the dark spots are artifacts due to contamination on the CCD). Note that 
some of the loops appear to cross each other. The number of loop crossing 
is much less than one would expect if the corona contained braided magnetic 
fields on scales £ > 5 arcsec. 

shows active region AR 1 1041 observed on 2010 January 25, 
starting at 11:02 UT. The observation used the Ti-poly filter, 
which is sensitive to plasma with temperatures greater than 1 
MK. In order to improve the signal-to-noise ratio, we added 
12 exposures with exposure time 0.5 1 s, taken over a period of 
50 minutes. Figure |2j5 shows a spatially filtered image where 
the high spatial frequencies have been enhanced to bring out 
the loop structures. 

The structure of the emission is complex, and multiple 
structures are superposed along the line of sight. There are 
some distinct coronal loops, but there is also a more diffuse 
emission component. The loop widths vary from about 3 arc- 
sec for the narrowest loops to about 20 arcsec for the widest 
ones. The wider loops show some variation in X-ray bright- 
ness across the loop on a scale of a few arcseconds, and these 
fine-scale structures appear to be aligned with the loop axis. 
However, our ability to observe these features is limited by 



the spatial resolution of XRT and (for the outer parts of the 
AR) by photon noise. 

The wrapping of bright loops around each other has occa- 
sionally been observed with the Transiti on Rej^ion and Coro - 
nal Explorer (TRACE) (see examples in Schri jver et al.ll 999^. 
but these cases are often ambiguous and most loops observed 
with TRACE do not show evidence for magnetic braiding. To 
search for braided fields with XRT, we looked for places in 
Figure |2] where two loops appear to cross each other There 
are only a few such sites within the observed active region. If 
the magnetic field were braided on observable scales (greater 
than 5 arcsec), one would expect to find many more such 
loop crossings. The few examples that can be found seem 
to involve loops that are well separated in height, and do not 
appear to be due to braided structures. For the wider coro- 
nal loops, it appears that the different threads within the loop 
are co-aligned to within a few degr ees, not 20° as predicted 
by Parker's original braiding model (lParkerll983t iPriest et al] 
l2002HKhmc huk 2006). We conclude there is no evidence for 
the existence of strongly braided magnetic fields in the corona 
on spatial scales of a few arcseconds or larger If there is 
braiding on the Sun, it must occur on small transverse scales 
(less than 5 arcsec) and involve small mis-alignment angles 
(at most a few degrees). 

2.2. The Coronal Heating Rate 

An important constraint on any model for coronal heating 
is that it must explain the average heating rate. The observed 
radiative and conductive losses of active regions imply a non- 
radiative energy flux into t he corona Fmech ^10^ erg cm"^ s"' 
jWithbroe & Novell 1 19771) . Assuming this energy enters a 
coronal loop at both ends and is distributed uniformly over the 
full length Lcor of the loop, the average volumetric heating rate 
Qcor = 2F,„gch/Lcor- For a loop with constant cross section and 
length Lcor ~ 50 Mm we obtain Qcor ^ 4 x lO"-' erg cm"^ s"' . 
Another method for estimating the heating rate is to use the 
scaling laws first developed by iRosner et al.l (1 19781 hereafter 
RTV): 



l.4xlO\pcorLcor/2) 



1.9 X 10" pI!^. 



1/3 



1/3 



^50 Mm 

Qcor~9.8xlO'plif{Lcor/2r^''' 



K, 



-5/6 



(1) 



= 1 .44 X 10"^ plif. ( ) erg cm " s 



50 Mm 



3-' (2) 



where T,„a^ is the peak temperature (in K), pcor is the coro- 
nal plasma pressure (in dyne cm"^), and Lcor is the loop 
length (in cm or Mm). X-ray observations indicate that the 
loops in the core of an active region have high temperature 
and pressure, Tmax ^ 2.5 MK and Pcor 2 dyne cm~^ (e.g., 
Saba & Str ong 119911; Brosius et al.l 119971; lyiiebarger et al.l 
200& ,Warren et al.ll2008i: ,Ugarte-Urra et al.1 12009'). Assum- 
ing a loop length Lcor ~ 50 Mm, the required heating rate is 
about 2.9 X lO"-' erg cm"-' s"'. 

Can the quasi-static braiding models explai n such heat- 
ing rates? In Parker's orig inal model (Parkerl lT97l IT983I) 
it is assumed that there exist distinct "flux tubes" that can 
be traced from the corona into the photosphere. These flux 
tubes are assumed to retain their identity for about 1 day 
as their footpoints are moved around on the photosphere, 
and thin current sheets develop at the interfaces between 
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the flux tubes (lBergeiill991lfT99l lWilm ot-Siiuth et al.l l2009t 



iBerger & Asgari-Targhill2009l) . To explain the observed heat- 
ing rates, large departures from the potential field must de- 
velop. Specifically, the angle 6 between the brai ded field B(r) 
and die pot ential field Bo(r) must be about 20° (Parker 1983, 
iPriest et aH i2002: KUmchuk .2006 ). otherwise the Poynting 
flux into the corona is insufficient to balance the observed ra- 
diative and conductive losses. These large deviations from the 
potential field are predicted to occur on transverse scales of a 
few megameters or larger, and therefore should be readily ob- 
servable in the fine structures of coronal loops. As discussed 
in Section im there is no observational evidence for braided 
fields on such scales. Therefore, it appears that Parker's orig- 
inal version of the braiding model (with relatively long-lived 
flux tubes braided on observable scales) is incompatible with 
coronal observations. 

In other braiding models the boundary flows are incom- 
pressible and vary randomly with time. This implies that 
any two p oints on the boundary te nd to separate from each 
other (e.g.. Ivan Balle gooiien "1988V hence there are no well- 
defined "flux tubes" that retain their identity for many hours. 
This type of model is consistent with the observation that 

S ihotospheric flux ele ments frequently split up and merge 
Berger & Titld Il996h . The magnetic fields produced in 
such models are more fragmented than those predicted by 
Parker's original model. Specifically, current sheets can de- 
velop anywhere within the volume, and the strongest cur- 
rents will develop in locations where the footpoint motions 
have the largest cumulative shear. Analytic studies predict a 
rapid "cascade" of magnetic energy to small transverse scales 
(|yan Ballegooiien 1985, 1986). In fact, the cascade proceeds 
so rapidly that the energy is dissipated before strong depar- 
tures from the potential field can develop. As a result, these 
simple cascade models have difficulty explaining the observed 
rate of coronal heating in active regions. The model predicts 
that the coronal heating rate per unit volume is given by 

Bl 2ulTo\nR 



(3) 



87r 3L2V2^ 

where Bo is the coronal field strength, mq is the rms veloc- 
ity of the footpoint motions, tq is the corre lation time, and 
j?m is the magnetic Reynolds number (see van Balle gooiienI 
Il986|) . For motions induced by the solar granulation, tq 5 
minutes, which is slow compared to the coronal Alfven travel 
time (ro ^ LjvA, where va is the Alfven speed, ^ 1000 
km s"'). Therefore, the magnetic field is expected to evolve 
through a series of nearly force-free equilibrium states. Note 
that Qcor is proportional to the product u^t^, which is essen- 
tially the photospheric diffusion constant, D = jUqToVTtt in 
this model. Therefore, measurements of D can be used as 
a constraint on the braiding model. Using D ss 250 km^ s"' 
(fPeVore et al.H l985). it was found that the heating rate pre- 
dicted by the cascade model is a factor 10 to 40 smaller than 
the heating rate observed in active regions. 

Numerical simulations of magnetic braiding driven by slow, 
random, incomp ressible foot p oint m o tions (e.g., 'Mikic et al. 
1989; Longcope & Strauss 1994; Hendrix & van Hoven 
"996: Hen drix et al. 1996; Galsgaard & Nordlund .1996.; 



Craig & Sn eyd 2005) predict heating rates that are basically 
consistent with equation (O. Therefore, like the above 
cascade model, these numerical models also cannot explain 
the observed heating rate when realistic values for the random 
walk diffusion constant D are inserted into the model. We 
conclude that neither Parker's original version of the braiding 



model nor its later modifications can adequately explain the 
structure and heating of coronal loops in active regions. This 
leads us to consider other types of footpoint motion, such 
as small-scale random motions inside the photospheric flux 
elements (see Figure[TJ. 

3. ALFVEN WAVE TURBULENCE MODEL 

A model for the dynamics of plasma and magnetic field in- 
side a coronal loop is developed. We consider only a thin 
tube of magnetic flux, corresponding to a single kilogauss 
flux tube in the photosphere. The tube extends from the pho- 
tosphere at one end of the loop, through the chromosphere 
into the corona, and back down into the photosphere at the 
other end. In the photosphere and chromosphere the flux tube 
is assumed to be vertically oriented. We assume that a sin- 
gle flux tube at one end of the loop is connected to a single 
flux tube at the other end, i.e., we ignore the fact that on the 
real Sun the photospheric flux concentrations at the two ends 
are uncorrelated and do not perfectly match up. Also, on the 
real Sun magnetic flu x elements frequently split up and merge 
with the ir neigh bors (iBerger et al.lll995i iBerger & Titielll996l; 
Ivan Ba llegooiie n~et alj Il998l) ! Such processes may lead to 
magnetic recon nection and may be important for coronal heat- 
ing (e.g.. .Furusawa & Sak ai 2000; Sakai & Fumsawa 2002). 
However, in the present model we neglect such effects, and 
we assume that the flux tube retains its identity for the dura- 
tion of the simulation. 

The tube has a length L, and the overall curvature of the 
tube is neglected. We use a coordinate system {x,y,z), where 
z is the coordinate along the tube axis (0 < z < L), and x and 
y are perpendicular to the axis. Note that near the "left" end 
of the flux tube (z « 0) the height in the lower atmosphere is 
given by z and gravity acts in the -z direction, but near the 
"right" end (z ~ L) the height is given by L-z and gravity is 
in the +z direction. Despite these differences, we will some- 
times refer to z simply as the "height" in the flux tube. The 
tube is assumed to have a circular cross-section with radius 
R{z), which is much smaller than its length L. To simulate the 
effects of solar convection interacting with the flux tube, we 
impose random footpoint motions on the field lines at the base 
of the photosphere (z = and z = L). These motions produce 
Alfven waves that travel along the flux tube, reflect due to gra- 
dients in Alfven speed, and generate turbulence via nonlinear 
wave-wave interactions. The dynamics of the waves inside 
the tube are described by the equation of motion. 



d\ 1 
p— =-Vp + pg+— (VxB)xB + D,, 

at 4-n 

and the magnetic induction equation, 
dB 



dt 



= Vx(vxB) + D„ 



(4) 



(5) 



Here p(r,f) is the plasma density, p{r,t) is the pressure, v(r,f) 
is the velocity, B(r, f ) is the magnetic field, g = gQ{z)i is the ac- 
celeration of gravity, and D,, and D„, are viscous and resistive 
dissipation terms. The magnetic field satisfies the solenoidal 
condition, V • B = 0. The velocity v can be split into parallel 
and perpendicular components: 



V = v_L + viiB, 



(6) 



where B(r, f ) is the unit vector along the perturbed magnetic 
field, and • B = 0. Taking the inner product of equation (|4) 
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with B(r, f ), we obtain 



dt 



dB 

~di 



+ B-(-V/7 + pg + D,). 



(7) 



The first term on the right hand side is the centrifugal force 
due to changes in the shapes of the magnetic field lines, which 
is an important driver of field-aligned flows. However, in the 
present paper we will neglect parallel flows and focus on the 
perpendicular motions of the plasma. Note that, according to 
equation (|5]l, the parallel velocity V[| has no direct effect on 
the evolution of the magnetic field. 

In this paper we use the so -called reduced MHD (RMHD) 
approximation dStraussll 19761) . which assumes that the mag- 
netic fluctuations are small compared to the background field. 
In subsection l3.1l we present a derivation of the RMHD equa- 
tions in the context of the present model. In the following 
subsections we describe the structure of the background at- 
mosphere, the imposed footpoint motions, and the numerical 
techniques for solving the RMHD equations. Subsection [33] 
describes some of the limitations of the present model. 

3.1. Reduced MHD Equations for a Non-Uniform Medium 
With Gravity 

In this section we derive the RMHD equations describing 
the nonlinear dynamics of Alfven waves in a flux tube with 
non-constant cross-section and with density varying by orders 
of magnitude alo ng the flux tube. The RMHD e quatio ns were 
first d erived by iKadomtsev & Pogutsd d 19741) and IStrauss 
dl976) for a uniform background field (also see Montgomery 
[1982; Hazeltine 1983), and the relationships between com- 
pressible MHD, incompressible MHD and reduce d MHD 
were extensively discusse d by IZank & MatthaeusI (Il992h . 
ISchekochihin et aP (l2009l) considered t he extension of the 
RMHD equations into the kinetic regime. [Bhattachariee et al.l 
([J998) included the effects of plasma pressure, spatial inho- 
mogeneities and parallel flows into the formalism, and they 
derived a more general set of four coupled equations. How- 
ever, they did not include the effects of gravity. In the ab- 
sence of gravity, the magneto-static equilibrium equation im- 
plies that the gradient of the plasma pressure is perpendic- 
ular to the mean magnetic field, Vpo -L Bq. In contrast, in 
the present case gravity plays an important role in the stratifi- 
cation of the plasma pressure poiz) and density poiz) in the 
photosphere and chromosphere, and at the axis of the flux 
tube Vpo II Bp. The r efore , the four-field equations derived by 
iBhattacharjee et aLl ([1998) cannot be directly applied to the 
present case. Also, many analyses of the RMHD equations 
use normalized dynamical variables, which makes sense when 
the background density is roughly constant, but not when po 
varies by several orders of magnitude within the system. This 
makes it necessary to discuss in some detail the assumptions 
underlying the equations used in the present work. 

Our starting point is the equation of motion (|4). The mag- 
netic and velocity fields are written as sums over mean and 
fluctuating components, B = Bo + Bi + • • • and = \±fl + 

Vx,i H , and similar for the parallel velocity, density and 

pressure. The background field Bo is non-uniform and varies 
on a spatial scale Hb, which is defined by 



Hb — Bq ( Bq • VBq 



(8) 



where Bo(r) is the background field strength, and Bo(r) is the 
unit vector along the background field. We assume that the 



background atmosphere is in static equilibrium (v^.o = and 
V|| = 0), and that the interior of the flux tube is current-free, 
V X Bo = 0, i.e., all currents are located at the interface be- 
tween the flux tube and its surroundings. Then equation ([4) 
yields Vpo = Pog, where pq{z) and poiz) are the unperturbed 
pressure and density as functions of height z inside the flux 
tube. 

We now assume that the radius Riz) of the flux tube is small 
compared to the length scale |//b(z)| of the magnetic field in- 
side the tube. Then we can define a smaU parameter, 

Riz) 



\Hb{z)\ 



«1, 



(9) 



and the background field Bo(r) can be approximated as 

Bo=Booz-5^(^+3'y) + 0(B()oeo), (10) 

where Boo(z) is the field strength on the tube axis {x = y = 0). 
The unit vector along the background field is given by 



1 9 

Bo = ^ " ^ (-*^ +3'y) + ^(^o) : 

ZHb 



(11) 



where we used Hb{z) ~ -Boo / {dB^) / dz). Flux conservation im- 
plies that B\x)R^ ~ constant along the flux tube. Note that 
V • Bo = as required, and that the unit vector varies over the 
cross-section of the flux tube: 

^-^=-^ + 0{el/R\ ^=-J- + 0{el/R). (12) 
dx 2Hb dy 2Hb ^ o/ ' y ' 

The above "thin tube approximation" has been used by many 
authors to study waves and instabilities in flux tubes (e.g., 
[Defouw 1976; Roberts & Webb 1978; Spruit 1981). 

We now consider the perturbations due to Alfven waves that 
are launched at the base of the photosphere (z = and z = L) 
and are reflected in the chromosphere and at the TR. The ve- 
locity amplitude u±{z) of the waves is assumed to be small 
compared to the Alfven speed va{z), so that the magnetic per- 
turbation Bi(r,f) is small compared to the background field 
Bo()(z)- Furthermore, we assume that the transverse length 
scale i± of the waves is less than the tube radius (£± < R) 
and is small compared to the parallel scale The latter is 
defined by 

i\\=mm(vAT,\HBlL), (13) 

where t is the typical time scale of the magnetic fluctuations 
(e.g., the Alfven wave period), \Hb\ is the length scale of the 
background field, and L is the loop length. Then we can define 
a second small parameter: 



e = max( — ,-— ) < 1, 



(14) 



and we can expand the magnetic and velocity fields in powers 
of e: 



B = Bo + Bi + 0(Booe^), 
vj. = v±,i + v±,2 + O(vAe^), 



(15) 
(16) 



where in general |Bi | = 0{Boo£), \y±.\ \ = O(vAe) and |v^.2| = 
0(vAe^). We assume that the main driver of parallel flows is 
the centrifugal force given by the first term in equation ([7). 
Then the parallel velocity V|| = ©(v^e^). 

In some derivations of the RMHD equations (e.g.. StraussI 
[19761 [Mont gomerv,.1982) it is assumed that the plasma flow 
is incompressible, V • v « 0, but others have argued that this 
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assumption cannot be take n for granted when the plasma 
beta is of order unity (e.g., IZank & Matthaeuslll992l) . Here 
(3 = Sttpo/Bo is the ratio of gas pressure and magnetic pres- 
sure. For the models presented in this paper, (3 ^ \ m the 
photospheric part of the flux tube (see section ITZt . It is not 
clear from the literature on the RMHD equations whether the 
assumption of incompressibility is still valid when (3^1 and 
the background medium is non-uniform. Therefore, we must 
first estimate the magnitude of V • v for our particular model 
and determine whether the assumption of incompressibility is 
still valid. Using equation the velocity divergence can be 
written as 

V-v= V-v_L + B-v(-^) . (17) 

The second term is of the order of vae^/^j., but the magnitude 
of the first term is unclear. To estimate this term, we insert 
expressions (fTSI l and (fT6l l into the momentum equation (|4]l, 
omitting the dissipative term, and then take the divergence, 
which yields 



(18) 



It follows that t he perturbation of the total pressure is of sec- 
ond order in e (IZank & Matthaeuslll992t iBhattachariee et al] 
fT998.) : 



Pi + - 



Bn • Bi 



47r 



= o 



n2 
47r 



(19) 



Also, taking the inner product of equation (|5]l with B/47r, we 
obtain 



1 

47r 



(v_L X B) • (V X B) 



-V I • I po— + '\7pi- Pig 



(20) 



where we used equation (|4|l, again without the dissipative 
term. Inserting expressions (ITSt . (IT6^ and (|T9t . we find 



dpi „ Bl ^ 



■ O 



47r^ I 



(21) 



Finally, the pressure is assumed to evolve adiabatically, 
dp/dt = -7pV • V, where 7 is the ratio of specific heats. In 
lowest order of e this equation becomes 

^ + y±-^pi+jpo\7-y^ = o(^po^e'y (22) 

where we used v • Vpo ^ VAe^\dpo/dz\- Adding equations 
(I2TI 1 and (I22I 1. and dividing by the factor (7/?o + Bo/47r), we 
find 



V-vi 



(23) 



i.e., the magnitude of V • is of third order in e. Inserting 
this result into ( |22] |. we find p[ = Oipoe^), so the pressure vari- 
ations are really of second order in e, and we can set = 0. 
Inserting pi = into equation (fT9] l. we find Bo -Bi = ©(BqoE^), 
so the first-order perturbation of the magnetic field is perpen- 
dicular to the background field, Bi _L Bq. Therefore, the ef- 
fects of compressibility can be neglected in the present mod- 
els, V • = 0, even though /3 ^ 1 in some regions of the 
model. This is due to the fact that the transverse motions of 



the waves are nearly horizontal and therefore along the planes 
of constant pressure poiz) inside the flux tube. 
The induction equation (|5]l can also be written as 



(9A 

— = VixB + V0, 
at 



(24) 



where A(r, t) is the vector potential (B = V x A), (/)(r, f ) is 
the electric potential, and we omit the dissipative term. From 
the above analysis it is clear that the Lorentz force is of or- 
der e^, so the componen t of electric cu rrent perpendicular to 
Bo is also of this order (IStrausslll997l) . Therefore, the first- 
order perturbation of the vector potential must be parallel to 
the background field: 

A,(r,f) = /i(r,f)B„(r), (25) 

where h is the magnetic flux function {h ~ £^e). Using V x 
Bo = 0, we find Bi = V/i x Bo, which is perpendicular to Bq as 
required. Then the total magnetic field is 

B(r,f) = Bo + Vi/zxBo + 0(Bgoe')- (26) 
Here is defined in relation to the background field, V = 
V - Bo(Bo • V), whereas vj^ is defined to be perpendicular to 
the perturbed field B. Inserting expression ( |25] l into (l24l l and 
taking the inner product with B, we obtain 

^ = ^B.V0=i2Bo-(V0+Vi0xVi/!), (27) 
at ts^ Hq 

where we use B • Bq = Bfy Taking the cross product of 
with Bo, and using the condition • B = 0, we find 



Vi(r,f)=-2(Vi<^xBo + ^Bo) + 0(yAe'), 

On 



(28) 



where ^ = -V^0 ■W±h. The two terms with and ^ are 
of order v^e and v^e^, respectively. Taking the divergence of 
equation ( |28] l. we obtain 



V-Vx = (Vi(^xBo + ^Bo)-V( -J l+-jBo-V^+0( 



1 



1 



va 3 
1 
(29) 

Estimating the magnitudes of the various terms, and using the 
fact that VBo is nearly parallel to Bo, we find that all terms 
are of the order of VAe^ /£±, so equation ( l29b is consistent with 
equation ( l23T l. Neglecting terms of order in equation ( l28b . 
the perpendicular velocity can be further approximated as 

Vi = Vi/xBo + 0(yAe2), (30) 

where / = <j>{r,t)/Bo{r). In this approximation the compo- 
nent of vj^ parallel to Bo is neglected. The quantity f(r,t) 
can be interpreted as the velocity stream function. Inserting 
the expression (j> = fBo into equation dZTI l. we obtain for the 
magnetic induction equation: 

5 = B„-V/+^+Bo-(Vi/xVx/i). (31) 
at Hb 

We now consider the equation of motion dUi. As discussed 
above, we can neglect the fluctuations in gas pressure and den- 
sity (pi = pi= 0), and we also omit the viscous term. Using 
equations (fTST i and ( |26] |, the curl of the magnetic field can be 
approximated as 

V X B = Bo • V( V_l/z) - V±/! • VBo - [V • ( V±/j)]Bo 

+ 0(BooeV^i) 

= Ba [Bo • V(Vi/z) + (2//B)"'V_L/i + aBo 



+ 0(BooeV^±), 



(32) 
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where a is defined by 

a(r,f) = -Vi/!. (33) 

Inserting expression (l32t into equation (0), we find for the 
acceleration perpendicular to Bq: 

- + V.VV 



Bo-V(Vi/i) + (2i/B)-iVi/! 



X Bo + aV^/z 



(34) 



where \'a{z) = Boo/V^Trpo- We now consider the parallel 
component of the vorticity equation. Using equation (|30] l. the 
parallel vorticity can be approximated as 

-vi/, 



B 



t)-(Vxv); 

1 ..2\ 



(35) 



and using v • Vv = V(5V ) + (V x v) x v, the inertial term in 
equation ( |34] l becomes 

Bo • [V X (V • Vv)] « Bo • (Vic^ X V^/). (36) 

The cross-product on the right-hand side of equation 
yields 



Bo 



V X I Bo • V(Vi/i) + (2//B)-iV_L/? X Boj 



= Bo-Va, 
(37) 

where equation (fT2l i is used several times, and Hb eventually 
drops out of the expression. Therefore, we obtain the follow- 
ing scalar form of the vorticity equation: 

diij 



^+Bo-(ViWxVi/) = v2 



Bo-Va + Bo-(V_Lax V^/i) . 

(38) 

Here va depends on position z along the flux tube, and Bq 
depends on x and y as described by equation (fT2] i. The mag- 
netic field strength Bo is nearly constant over the cross-section 
of the tube and can be approximated by its value on axis, 
Bo ~ Boo(z)- Therefore, in the remainder of this paper we will 
write the magnetic field strength as Bo(z)- 

The last step is to replace with = V-z(z • V) in 
the definitions of a and uj, and in the nonlinear terms of 
equations ( l3Tl i and (|38] |. but not in the linear terms where 
quantities are differentiated along the background field. Then 
V^i « /dx^+d^ /dy^, and the dynamical equations (l3Tl i and 
(|38] | can be written as 

^=Bo-V/+-^ + [/,/;], (39) 
of Hb 

|^=-[c^,/] + vi{Bo-Va+[«,/i]}, (40) 

where [•••,•••] is the bracket operator: 

da db da db 
dx dy dy dx 

We also derive an alternative form of the dynamical equations. 
Taking the second derivative of equation (|39] l, we find 

da 



(41) 



dt 



-Vi(Bo.V/+[/,/i]) + 



Hn 



= Bo-Vt^+[w,/!] + [/,a]-2 



'df dh 


-2 


'df dh 


dx ' dx 




dy ' dy 



(42) 



where we used equation ( fT2] l to evaluate the horizontal deriva- 
tives of Bo (note that Hb is eliminated from the equation). In 
analogy with the Elsasser variables (Elsasser 1950), we define 

f±=f± vaIi, and w± ee w ± v^a, (43) 



and inverting these relations, we can write f = (f+ + /_) /2, 
h = (/+ - /-)/(2v^), uj = (a;+ + aj-)/2 and a = {uj+-uj-)/(2va). 
Therefore, combining equations (l40t and (l42t . we obtain the 
following equations for uj±: 



duj± 
"df 



= ±va Bo • Va;± ■ 



dVA 

-va— — a- 
dz 



[oj±,f^]±M, (44) 



where the nonlinear term J\f is defined by 



'dU df. 




'dU dfT 


dx ' dx 


+ 


dy ' dy 



(45) 



The equations for uj+ and u- describe inward and outward 
propagating Alfven waves, respectively (T^o directions); the 
term with dvA/dz describes linear coupling between these 
waves; and the last two terms describe nonlinear coupling. 
In Appendix A we demonstrate that equation (l44l i is consis- 
tent with earlier formulations of the wave transport equations 
based on the Elsasser variables. The numerical methods used 
in solving these equations are discussion in section \3A\ and 
Appendix B. 

3.2. Background Atmosphere 

We now describe the undisturbed state of the flux tube be- 
fore any waves are injected. In the closed field case there are 
two "lower atmospheres" located at each end of the flux tube. 
The chromosphere-corona TRs are located at positions z = Ztr 
and z = Lcor+ZTR, where ztr is the height of the first TR, and 
Lcor is the coronal loop length. The total length of the tube 
is L = Lcor + '^ZTR- The gas pressure pq(z), density poiz) and 
magnetic field strength Bq{z) are functions of z only, i.e., we 
neglect variations of these quantities in planes perpendicular 
to the flux tube axis. Furthermore, the functions po(z), poiz), 
R(z) and Bo(z) are assumed to be symmetric with respect to 
the mid-point of the loop (z = L/2), so the two halfs of the 
loop are identical. 

The structure of the lower atmosphere is based on a model 
of a facular region (i.e., very bright plage region) called model 
P, which was first developed by Fontenla et al. (1999) and 
more recently discussed by Fo ntenla et al.l (120061 l2009). The 
temperature 7b(z) and molecular weight ^q(z) are specified as 
function of height to be in rough agreement with this model. 
At the base of the photosphere, we assume Bo(Q) = Bphot = 
1400 [G] and po(0) = 3 x 10"* dyne cm"^, which is the internal 
pressure of the flux tube; then the total pressure ipQ + B^/diTT) 
is about 1.1 X 10^ dyne cm"^, consistent with model R The 
gas pressure pq{z) and density pdz) as functions of height 
are determined by solving a modified hydrostatic equilibrium 
equation, dp^ jdz = -pogeff, where ggfj is the gravitational ac- 
celeration corrected for the effects of turbulent motions. The 
pressure at the TR is given by 



Po(ztr)= 1.833exp 



Ztr-LS 
" 0.3278 



dyne cm 



(46) 



where ztr is in Mm. The magnetic field strength in the lower 
atmosphere is given by 



Boiz) = 



(B 



phot 



Poiz)- 



Po(ztr) 2 



Po(0)-po(ztr) 



1/2 



forO<z<Z7 



(47) 

. is the field strength at the coro- 
nal base (i.e., at z = ztr)- In the photosphere po ^ po{ztr) and 
Bo ^ Bcor, so that the ratio of gas- and magnetic pressures is 



where Bp^ot = 1400 G, and B^ 
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approximately constant: 



87rpo(z) 
Bliz) 



■QA. 



(48) 



At larger heights in the chromosphere Bq{z) approaches the 
coronal value, and /3(z) decreases with height to values well 
below unity. 

In the TR the temperature Tq{z) increases rapidly with 
height from about 10"* K in the upper chromosphere to about 
10^ K in the corona. We do not attempt to resolve this tem- 
perature structure, and treat it as a discontinuity. For closed 
loops we use the following profile for the temperature in the 
corona: 

mz) = T,na.[l-0.^u\z)f'\ (49) 

where u(z) = -I + 2(z - Ztr) / Lc„r, which lies in the range 
-1 < M < +1, and T„u,x is the peak temperature in the loop 
as predicted by the RTV scaling law, equation ([U. It fol- 
lows that To(ztr) = 0.63 1 Tmax, so the temperature in the upper 
TR is a fixed fraction of the peak temperature in the corona. 
The above profile is similar to that predicted for a conduction- 
dominate d loop with constant cross-section an d unifor m heat- 
ing (e .g., IVeseckvetanil979l: iMartens etar.i.2000; MartensI 
I2OIOI) . Alth ough it is not dif ficult to include the effects of 
gravity (e.g.. lSerio et al.lfr981h . in this paper we neglect grav- 
ity in the corona, so that the gas pressure poiz) is constant 
along the coronal part of the loop. The pressure is continuous 
across the TR, so the coronal pressure pcor = Po{ztr)- There- 
fore, the coronal pressure is determined by the height ztr of 
the TR. The plasma density in the corona is given by 



Poiz) = po(ztr) 



1-0.8m^(z) 
0.2 



-2/7 



(50) 



where po(ztr) is the density at the coronal base, which is com- 
puted from Per and ToizrR)- For closed loops the coronal field 
is approximated as 

Boiz) = B„„. { 1 + (r- m-u\z)]y' , (51) 

where F is the areal expansion factor; for open fields we use 
Boiz) = Bcor = constant in the corona. In either case the mag- 
netic flux $ (= ttR^Bq) is constant along the flux tube, so the 
tube radius is given by 



R(z) = Rpho,jBpho,/Bo{z) 



(52) 



Here Rphot = 100 km is the tube radius at the base of the pho- 
tosphere (z = and z = L). The radius R in the corona depends 
on the parameters Bcr and F, and is different for the different 
models. 

We consider a reference model with coronal field strength 
Bcor = 50 G, expansion factor F = 1, and TR height ztr =1-8 
Mm, which corresponds to a coronal pressure pc„r = 1.833 
dyne cm"^. The coronal loop length Lc„r = 49.6 Mm, which 
is the closest one can get to 50 Mm with the chosen time step 
Afo = 0.746 s (see Appendix B for details). The total loop 
length is L = 53 .6 Mm. Figure[3]shows various quantities plot- 
ted as function of position along the flux tube for this model. 
Positions are given in terms of the Alfven wave travel time 
from the left footpoint (z = 0): 



T(Z): 



dz 



(53) 



/o va{z) 

Figure [3^ shows the relationship between z and r for the 
reference model; the photospheric footpoints are located at 



r(0) = and r(L) = 190.9 s, and the corona is located in the 
region 81.3 < r < 109.6 s. The other panels in Figure|3]show 
the temperature Tq, density po, magnetic field strength Bo, flux 
tube radius R, the absolute value of the magnetic scale height 
Hb (defined in equation (O), and the Alfven speed as func- 
tions of Alfven travel time r from the left footpoint. Note that 
when expressed in terms of r, the corona is only a small part 
of the computational domain. 

Figure[3^ shows that \Hb\ < in two intervals, 30 < r < 48 
[s] and 140 < r < 158 [s], which correspond to the tem- 
perature minimum regions at the two ends of the loop. In 
these regions the quantity eo defined in equation ^ is greater 
than unity, so the thin tube approximation is no longer valid. 
Clearly, the thin tube and reduced MHD approximations dis- 
cussed in section [TTI provide only a crude description of the 
magnetic structure and wave dynamics in the lower atmo- 
sphere. A proper treatment will require full MHD simula- 
tions, and is beyond the scope of the present work. 

The present model predicts braiding of the coronal field 
lines on a transverse scale less than the tube radius, Rcor ~ 529 
km. This radius is less than the resolution limits of present- 
day X-ray telescopes. Therefore, we should expect that the 
predicted coronal structures will be difficult to observe. 

3.3. Photospheric Footpoint Motions 

The Alfven waves are produced by footpoint motions im- 
posed at the two ends of the flux tube, z = and z = L. In 
reality these motions may distort the shape of the flux tube as 
indicated in Figure [T] but here we use a simpler approach in 
which the motions are assumed to be confined to a circular 
area x^ + y^ < R^pij„, at z = 0, and similar at z = i. The veloc- 
ity \(x,y,Q,t) at these boundaries can be written in terms of 
polar coordinates (r, (p), where r is the distance from the flux 
tube axis and ip is the azimuth angle in the (x,y)-plane. We 
assume that the radial component of velocity vanishes at the 
tube wall, VriRpiwi , V, 0, = 0, so that the circular shape of the 
cross-section is preserved. We also assume that the motions 
are horizontal and incompressible: 



v(r,(^,0,f) = V/xz, 



(54) 



where f{r,Lp,Q,t) is the stream function at z = 0, and sim- 
ilar at z = i. As described in Appendix B, functions on a 
circular area can be decomposed into orthogonal basis func- 
tions Fi(£^,(p), where ^ is the relative distance from the tube 
axis (^ = r/R) and index / enumerates the basis functions 
(/ = 1 , • • • , A^). We use a relatively small number of basis func- 
tions (A^ = 92); the functions are shown in Figure [13] 

In this paper we assume that the footpoint motions have a 
pattern consisting of two counter-rotating cells. This pattern 
can be described as a superposition of two modes with az- 
imuthal mode number m= I. For the particular set of basis 
functions used in the present paper, these driver modes have 
indices / = 7 and / = 8 (see Appendix B for details), and are 
shown in the top row of Figure [T3] (seventh and eighth image 
from the left). Both modes have the same radial dependence 
given by jQ{a±r/Rpho,), where J()(x) is the zeroth order Bessel 
function of the first kind, and a± is the dimensionless perpen- 
dicular wavenumber, which equals 3.832 for these particular 
modes (the first zero of the Bessel function). However, the 
azimuth dependence of the two modes is different: the mode 
with / = 7 is proportional to cos (p, while the one with ; = 8 is 
proportional to sin (p. The imposed stream function at z = 
can then be written as a superposition of the two modes: 

f(r,^,0,t) = Mt)Fj(r/Rpho„p) + Mt)Fsir/Rpho„^). (55) 
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Fig. 3. — Reference model for a coronal loop and the two "lower atmospheres" at the two ends of the loop. Vaiious quantities are plotted as function of 
the Alfven wave travel r to a given point along the loop: (a) position z(r) along the loop as measured from the left footpoint; (b) temperature Tq; (c) mass 
density pq; (d) magnetic field strength Bq; (e) flux tube radius R (full curve) and magnetic scale height \Hb\ (dashed curve); (f) Alfven speed v^. The two 
chromosphere-corona TRs are located at r = 81.3 [s] and r = 109.6 [s]. 



Here fi(t) and /8(f) are the mode amplitudes, which vary ran- 
domly with time t and are not correlated with each other The 
vertical component of vorticity at z = is then given by 

u(r, (^,0,0 = -Vi/ = (3.S32/R„ho,ff(r, 0, f )• (56) 

The above time-dependent pattern simulates the intermixing 
of the plasma within the flux tube due to motions imposed by 
the surrounding convective flows. 

The random variables fj(t) and /§(?) are constructed as fol- 
lows. For each variable, we first create a normally distributed 
random sequence /(f) on a grid of times covering the entire 
length of the simulation (f,„at = 3000 s). Then the sequence 
is filtered in the Fourier domain using a Gaussian function, 
G(P) = exp[-(ro!>)^], where P is the temporal frequency (in 
Hz) and tq is the specified correlation time (for the reference 
model To = 60 s). The filtered sequence is renormalized such 
that the rms vorticity of each mode equals a specified value, 
ojQ. The root-mean-square (rms) velocity of the combined 
pattern of the two modes is Avrms = ViR phoii^o/ 3. S32. The 
quantities tq and are free parameters of the model. 

3.4. Numerical Methods 

The techniques used for solving the RMHD equations are 
described in Appendix B, and the energy equation is discussed 
in Appendix C. The transverse dependence of the waves is de- 
scribed using a spectral method, i.e., all scalar functions are 
written as sums over 92 discrete modes. The nonlinear terms 
in the equations are represented by a matrix Mty,- that cou- 
ples the different modes, and cause transfer of energy from 
low to high wavenumber The modes with the highest trans- 
verse wavenumbers are artificially damped, which describes 
viscous and resistive processes on small spatial scales. The 



damping rate is given by equation (IB 13b with I'max = 0.7 s~'. 
For the z-dependence of the waves we use finite differences; 
for example, in the model shown in Figure |3] there are 259 
points along the tube. To accurately simulate the wave prop- 
agation, we use a grid that is uniform in Alfven wave travel 
time t(z) with a grid spacing At equal to the time step Afo 
of the simulation (for the reference model, Afo = 0.746 s). At 
the chromosphere-corona TR the density pq(z) is discontinu- 
ous, which is represented by two grid points at the same posi- 
tion ztr but with different densities. The discontinuity causes 
wave reflections that are described in terms of reflection and 
transmission coefficients (see Appendix B for details). 

The RMHD model is valid only when AB„„s/Bo < 1, 
where AZ?„„5 is the rms value of the transverse magnetic field 
fluctuation. We will find that this condition is only marginally 
satisfied. Therefore, the present model can describe only 
some of the dynamical processes that occur in the chromo- 
sphere and corona. In particular, the model does not describe 
field-aligned flows such as spicules. 

3.5. Limitations of the Model 

The present model for chromospheric and coronal heating 
has several drawbacks and limitations. The code uses an ex- 
plicit numerical scheme, which makes it difficult to simulate 
waves in models with high coronal field strength. On the real 
Sun the field strength in active region loops in the low corona 
is 100 - 500 G, but in the present paper we must limit our- 
selves to Bcor < 200 G. Also, the model includes only a lim- 
ited number of wave modes (see Figure [13), and the time step 
is relatively large (Afo > 0.1 s), so the turbulent spectrum is 
not well resolved. 

The model considers only a single magnetic flux tube, and 
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does not allow for splitting and merging of m agnetic flux el- 
ements at the two ends of the coronal loop (iBerger & Titlel 
|1996l). It is unclear to what extent the chromospheric and 
coronal heating rates and the spatial distribution of the heat- 
ing are affected by this approximation. Perhaps more impor- 
tantly, the present RMHD model does not allow for feedback 
of the heating on the background atmosphere. In the present 
version there are no field-aligned flows, and the temperature 
and density are assumed to be constant in time, so we do 
not simulate the dynamic response of the atmosphere to heat- 
ing events. This also means that we cannot make predictions 
regarding the temporal variability of EUV and X-ray emis- 
sions from the modeled coronal loops. In particular, we do 
not simulate spicules or similar dynamic events, even though 
such events may in fact be driven by nonlinear Alfven waves 
(see Matsumoto & Shibata 2010). The plasma density in our 
model is assumed to be constant in planes perpendicular to the 
tube axis, so there are no variations in density on scales less 
than the tube radius. Therefore, the present model has little to 
say about the multi-thermal structure of coronal loops. 

Finally, the model does not include the effects of phase mix- 
ing jH eyvaerts & Priest 1983) or resonant absorption (e.g., 
iDe Groof & Goossens,,2002.) . which depend on the variation 
of Alfven speed in the plane perpendicular to the tube axis. 
In the present model both Bo and po are constant in the {x,y) 
plane. 

4. RESULTS 
4.1. Reference Model 

We first construct a reference model with coronal field 
strength Bcor = 50 G, expansion factor F = 1, and coronal loop 
length of about 50 Mm. More precisely, the coronal loop 
length Lcor = 49.6 Mm, the transition-region height ztr =1.8 
Mm, and the total tube length L = 53.2 Mm. The footpoint 
motions have a correlation time tq = 60 s, each of the driver 
modes has a vorticity loo = 0.04 s"', and the rms velocity is 
Avrmi = 1 -48 km s"' . The latter is reasonable compared to the 
convective velocities of several km/s expected to exist in the 
downflow lanes below the photosphere. The TR height corre- 
sponds to a coronal pressure pcor = 1 -833 dyne cm~^, which is 
typical for hot loops found in active regions, and equation dU 
yields a peak temperature T^^ix = 2.3 MK. 

The vorticities of the driver modes are shown as function 
of time in Figures |4^ and ^ for the left (z = 0) and right 
(z = L) footpoints, respectively. These random footpoint mo- 
tions create Alfven waves that propagate upward along the 
flux tube. Initially, only the two driver modes are present. 
The decrease of density with height causes the velocity am- 
plitude of the waves to increase with height. In the chromo- 
sphere, part of the wave energy is reflected back down due 
to the increase of Alfven speed with height. Even stronger 
reflection occurs at the chromosphere-corona TR, where the 
Alfven speed suddenly increases by about a factor 15. These 
reflections produce a pattern of counter-propagating waves 
in the photosphere and chromosphere at the two ends of the 
loop. The wave amplitudes in the chromosphere soon build 
up to about 20 - 30 km/s, consistent with the observations by 
iDe Pontieu etalj (.2007b.) . 

Due to evolution of the driver modes and time de- 
lays in the reflection, the inward and outward propa- 
gating waves at a given height z have different spa- 
tial distributions in the {x,y) plane. Furthermore, 
these patterns are superpositions of multiple modes with 
different perpendicular wavenumbers. Such counter- 



propagating waves with different spatial patterns and multi- 
ple modes are subjec t to stron g non h near interactions (e.g., 
'Shebahn et al.' 1983"; 'Hi gdonI Il984t lO ughton & Matthaeusl 
1995; Goldreich & Sridh aHll995l 119971' Bhattachariee & 
2001; Cho et all 120021: lOughton et a l. 2004). In our RMHD 
model, these interactions are due to the bracket terms in equa- 
tions ( [39] l, (|40] | and ( |44] |. Basically, the inward propagat- 
ing waves are disto rted by the outward p ropagating waves, 
and vice versa (e.g. IChandran et ani2009h . These distortions 
are large because the fluid displacements are comparable to 
the transverse scale of the waves. For example, for chromo- 
spheric waves with velocity amplitude of 10 km/s acting over 
a period of 50 s, the transverse displacement is about 500 km, 
equal to the transverse scale of the driver modes. As a re- 
sult of these nonlinear interactions, other wave modes with 
smaller spatial scales are excited (see Figure [T3]for a display 
of the 92 wave modes used in the present model). After about 
200 s from the start of the simulation, a well-developed spec- 
trum of Alfven waves has formed and dissipation of the high- 
wavenumber waves becomes significant. This dissipation is 
due to a combination of viscous and resistive diffusion effects, 
so the model includes the effects of magnetic reconnection. 
Such turbulent dissipation of waves in the lower atmosphere 
continues throughout the simulation. 

A small fraction of the wave energy is transmitted through 
the TR into the corona. Energy is injected into the coronal 
loop at both ends, producing counter-propagating waves in 
the corona. As in the chromosphere, the counter-propagating 
waves significantly distort each other because (1) they have 
different spatial patterns and (2) the fluid displacements are 
comparable to the transverse scale of the waves. Therefore, 
the waves produce turbulence in the corona, and there is a 
continual cascade of energy to smaller transverse scales. As a 
result, part of the Alfven wave energy is dissipated in the coro- 
nal part of the loop, again due to combination of viscous and 
resistive effects. We find that the rms velocity amplitude of 
the waves in the corona is similar to that in the upper chromo- 
sphere, but the energy dissipation rate in the corona is much 
lower because of the lower coronal density. 

In the present model, all of the dissipation occurs 
via Alfven wave turbulence, both in the chromosphere 
and in the corona. Phase mixing and resonant ab- 
sorption of Alfve n waves (e.g. , iHeyvaerts & Priesd Il983t 
Poedts et a l.l Il989t lOfman et alJ fl994l: lErdelyi & GoossensI 



1 1 995£ iDe Groof & Goossensll2002i)" play no role because these 
processes require variations in Alfven speed across the field 
lines, which are not included in the present model. A key 
feature of the model is that the photospheric footpoint mo- 
tions include more than one dri ver mode, not just t h e torsi onal 
mode as in the simulations by lAntolin & Shibatal (1201 Ol) and 
iMatsumoto & Shibata (20fO). Specifically, we use two driver 
modes (modes k = l and ^ = 8 in Figure [T3T l with azimuthal 
mode number m= 1 . If the system is driven using only a sin- 
gle mode (e.g., k = 7), all nonlinear terms in equations (l40t 
and ( |39] | vanish. In this case all energy remains in the driver 
mode, and no turbulence develops, as we have verified in test 
calculations. Hence, for Alfven wave turbulence to develop it 
is essential that the footpoint motions include multiple driver 
modes that have nonlinear couplings with high wavenumber 
modes. This coupling can be indirect; for the chosen m=\ 
modes, the coupling runs via the m = modes. 

Figure|4j; shows the spatially averaged energy density of the 
waves E{t) in the corona as function of time (average between 
z = Ztr Sind z = ZTR+Lcor)- This quantity is the sum of magnetic 
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Fig. 4. — Vaiious quantities plotted as function of time for the reference model: (a) imposed voiticity at the left footpoint (z - 
k = l {dashed curve); (b) similar for the light footpoint = L); (c) energy density E in the corona; (d) coronal heating rate, Qcor- 



0) for k = 6 (full curve) and 



free energy E,„ag(t) and kinetic energy Ekmit), but the mag- 
netic energy dominates. Figure HJl shows the spatially aver- 
aged heating rate QcoAt) per unit volume in the corona. Note 
that both the energy density and the heating rate vary strongly 
with time. Significant heating starts only after about 200 sec- 
onds when the waves have become sufficiently turbulent. Dur- 
ing the latter part of the simulation, the average heating rate in 
the coronal part of the loop is 2.98 x lO""* erg cm~^ s"^. This 
is consistent with the second RTV scaling law, equation (|2|i, 
which yields Qcor = 2.9 x lO"-' erg cm"^ s"'. Therefore, the 
reference model describes a coronal loop in which the average 
rate of plasma heating is balanced by radiative and conductive 
losses. The model represents the kind of hot, dense loops typ- 
ically found in active regions. We conclude that such loops 
can be heated by Alfven wave turbulence generated inside the 
magnetic flux tubes by footpoint motions with rms velocity of 
about 1.5 km/s. 

Figure |5] shows various quantities as function of position 
along the flux tube, averaged over the cross-section of the flux 
tube {x and y) and over the time interval t = [800, 3000] s. Po- 
sitions are given in terms of the Alfven travel time r(z) in sec- 
onds. Figure|5t shows the kinetic and magnetic energy densi- 
ties, Eki„(z) and E„,ag{z), and their sum E(z). Figure|5}l shows 
similar plots for the kinetic and magnetic heating rates, Qkmiz) 
and Qmagiz), and their sum Q(z)- These quantities are dis- 
continuous at the TR. In the deep photosphere (t < 20 s and 



T > 170 s) and in the corona (8 1.3 < r < 109.6 s) the magnetic 
energy dominates, but in the chromosphere Eki„ > E,„ag', this is 
similar to what happens in open-field models with non-WKB 
wave re flection dominated by low-frequency waves (see Fig- 
ure 6 in ICranmer & van Ballegooijenll2005l) . Despite these 
differences in wave energy density, magnetic heating domi- 
nates over viscous heating almost everywhere in the model 

iQmag ^ Qkin)- 

The middle and right panels in Figure |5] show the rms val- 
ues of the velocity, vorticity, magnetic field fluctuation, and 
a parameter All four of these quantities are continuous at 
the TR. The rms velocity Av,,,,, and rms magnetic field fluc- 
tuation ABrms are defined in equations (IClb and (IC2b . and 
are further averaged over time. Note that the velocity and 
vorticity have their peak values at the mid-point of the loop 
in the corona, Avv/ 



37 km s and cj^ 



0.52 S-' at 



T = 95.4 s (see Figures ^ and |5};). Also note that the ve- 
locity at the mid-point is about 60% larger than that at the 
two TRs (vertical dashed lines in This is due to the 

fact that the waves in the corona are reflected back and forth 
between the TRs several times before they decay, creating 
standing waves with nodes near the TRs. The velocities of 
20 - 40 km/s found here for the corona are similar to the 
nonthermal velocities of 20 - 60 km/s found in observations 
of spectral line widths in active regions (e. g., iDere & MasonI 
119931: 1 Warren et all 120081: iLi & Dindl2009h and on the quiet 




Fig. 5. — Model for Alfven wave turbulence in a coronal loop. Various quantities are plotted as function of position along the flux tube for the reference model: 
(a) kinetic and magnetic energy densities, and their sum (b) rms velocity Avrms', (c) rms vorticity u],ms, (d) kinetic and magnetic heating rates, and their 
sum QizY, (e) rms magnetic fluctuation ABrms', (f) rms twist parameter arms- These quantities are averaged over the cross-section of the flux tube and over time. 
Positions along the loop are given in terms of the Alfven travel time t(z). The corona is located in the range 81.3 < t < 109.6 s, which is indicated by vertical 
dotted lines. 

erence model at time t = 2702 [s] in the lower atmosphere 
and in the corona, respectively (the vertical scales of these 
images are compressed by different factors). The field 
lines are traced from randomly selected points at height 
z = in the photosphere. The field lines are significantly 
distorted due to the Alfven waves that travel up and down 
the flux tube with a range of transverse wavenumbers. Two 
movie sequences of such images are available in the on-line 
version of the manuscript. These movies show the evolution 
of the magnetic field over a period of 298 seconds (from 
t = 2702 [s] Xo t = 3000 [s]), and are traced from footpoints 
that move with the flow. The coronal field lines are to 
some degree twisted and braided around each other (see 
Figure |7j)), but these structures are highly dynamic and 
change on a time scale of seconds. Therefore, the system 
is not in a force-free state, and is best described as Alfven 
wave turbulence. The effects of such turbulence on the 
solar corona ha v e been modeled previous l y for both open 
jHollwed 119861: IZhou & M atthaeusi 119901; 'Mat thaeus etal] 
1999t iDmitruk et all 12001; Dmitruk & MatthaeusI boW 



Sun (IChaeetal.lll"998t iMcIntosh et ani2008h . This suggests 
that the modeled waves have more or less the correct ampli- 
tude. The magnetic fluctuation ABrmAz) and twist parameter 
Oirmsiz) shown in Figures and |5f have their peaks in the 
lower atmosphere. This is due to the large gradients in Alfven 
speed in the chromosphere and TR, which cause strong wave 
reflection and produce a patterns of nearly standing waves in 
the lower atmosphere. The non-constancy of armAz) implies 
that the system is far from a force-free equilibrium state. The 
peak value of AB„„j/Z?() occurs in the low chromosphere and 
is about 0.3, so the RMHD approximation is only marginally 
satisfied. 

Figure |6] shows cross-sections of the tube at various posi- 
tions along the loop at the end of the simulation (f = 3000 s). 
The top row shows the velocity stream function f(x,y,z), the 
second row shows the vorticity uj(x,y,z), the third row shows 
the magnetic flux function h(x,y,z) and the bottom row shows 
the twist parameter a(x,y,z)- The different columns corre- 
sponds to different positions z along the loop, and are labeled 
with the Alfven travel time r(z) from the left footpoint. The 
two TRs are located at r = 81.3 s and r = 109.6 s. The upper 
left and upper right panels show the footpoint motions that 
drive the system. The vorticity uj{x,y,z) and twist a{x,y,z) 
exhibit small-scale structures that are produced by nonlinear 
interactions, as described above. A movie sequence of such 
images is available in the on-line version of the manuscript. 
This sequence covers the last 298 seconds of the simulation 
and shows that the system is highly turbulent. The waves in 
the corona have smaller spatial scales and evolve more rapidly 
than those in the lower atmosphere. 

Figures |7^ and [TJ? show magnetic field lines in the ref- 



Cranmer & van Ballesooiien 2003t i2005i: iCranmer et al.l 
2007; Verdini & Velli 2007) and clo sed magnetic fields 
(Hevvaerts & Priest 1984, 1992; Lonacope & Strausjl 
1994; Dmitruk & Gomez 1997; Buchlin & Velli 2007t 
Rapp azzo et al. 1 120081) . The present work demonstrates that 
Alfven wave turbulence can occur both in the chromosphere 
and in the corona, and can develop even when the photo- 
spheric footpoint motions occur on very small spatial scales 
(£± < 100 km). The model can quantitatively explain the 
observed heating rates in active regions. 
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Fig. 6. — Spatial distribution of various dynamical quantities in the reference model at time t = 3000 [s]: velocity stream function /, vorticity ui, magnetic 
flux function h, and twist parameter a. The different columns correspond to different positions along the flux tube, and are labeled with the Alfven travel time 
t{z). Each panel shows the normalized distribution of the relevant quantity as function of transverse coordinates .v and )' (see Figure |5]for information about 
normalization). A movie sequence is available in the on-line version of the manuscript. 





(a) 



Fig. 7. — Magnetic field hnes in the reference model at time t = 2701.7 [s], viewed from an angle of 30°. (a) The lower atmosphere up to the height of the 
first TR, ztr = 2 Mm. The starting points of the field lines are randomly distributed inside the flux tube at height z = (cylinder base). The radius of cylinder is 
0.53 Mm, and the vertical scale of the image is compressed by a factor 1 .9. (b) Continuation of the same field lines into the coronal part of the loop. The actual 
length of cylinder is 49.6 Mm, so the vertical scale of the image is compressed by a factor 47.3. Two movie sequences are available in the on-line version of the 
manuscript. 



4.2. Power Spectra 

The spatial power spectra for velocity and magnetic field 
fluctuations are defined in equation (|C4t . We computed 
such spectra for both the standard reference model (described 
above) and for a modified version in which the spatial res- 
olution of the model was slightly increased. Specifically, 
the maximum perpendicular wavenumber was increased 
from 20 to 26, which causes the number of modes to increase 
from 92 to 158. Also, the maximum damping rate was in- 
creased by a factor (26 /20)^, so that the damping rate vj; for 
the low wavenumber modes (a^t < 20) is the same as that in 
the standard model [see equation (IBOH . We found that the 
increased spatial resolution has little effect on the power in 
the low wavenumber modes, therefore, the standard model is 



adequate for most purposes. Nevertheless, in the following 
we show results from the modified version with a„ax = 26. 

Figures[8^ and[8]5 show velocity and magnetic power spec- 
tra binned in intervals of the dimensionless wavenumber a± 
for four heights in the reference model. Specifically, the ve- 
locity power in bin n is given by Py ,, = J^k^v.k, where Py^k is 
defined in equation (IC4b and the sum is taken over all modes 
with ak in the range nAa < ak < (n + l)Afl (n = 0, • • • , 12). 
Here Aa = 2 is the bin size in wavenumber space. A sim- 
ilar expression holds for the magnetic power spectrum Pb „. 
The results shown in Figure |8] were derived from the last 
800 time steps of the simulation (597 seconds). The solid 
curve in Figure |8^ shows the velocity power spectrum at the 
base of the photosphere (z = 0), and is dominated by the two 
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Fig. 8. — Power spectra and related quantities for the reference model: (a) velocity power spectra, (b) magnetic power spectra, (c) average frequency lu of 
the velocity fluctuations, (d) parameter describing the degree of nonlinearity of the waves (<^j_ = /rxi'x/i^). These quantities are plotted as function of the 
dimensionless perpendicular wavenumber ax - The different curves indicate different positions along the loop: photospheiic base (z = 0, solid curves, diamond), 
temperature minimum region (z = 0.50 Mm, dotted curves), chromosphere (z = 1.52 Mm, dashed curves), corona (z = 25.2 Mm, dash-dotted curves). 



driver modes with = 3.832. As we move to larger heights, 
the spectrum is broadened and the total power is increased. 
At height z = 1 .52 Mm in the chromosphere, the turbulence 
has generated a broad distribution of modes extending up to 
the maximum available wavenumber {a,nax = 26). The dash- 
dotted curve in Figure [8^ shows the power spectrum near the 
mid-point of the coronal loop (z = 25.2 Mm) where the level 
of turbulence is further enhanced. Figure |8j) shows the cor- 
responding curves for the magnetic power spectrum. At the 
base (z = 0, solid curve) the magnetic spectrum extends over 
a broad range of perpendicular wavenumbers and is very dif- 
ferent from the velocity spectrum at that height. The high 
wavenumber part of this spectrum is due to the downward 
propagation of waves produced in the chromosphere. The 
magnetic fluctuations in the corona are much smaller than 
those in the lower atmosphere at all wavenumbers. 

For each wave mode k and height z, we can also determine 
the temporal power spectrum Pk(uj,z) of the velocity fluctua- 
tions. Here ui is the wave frequency in radians per second. We 
compute this power spectrum by taking the Fourier Transform 
of the velocity stream function fii{z,t) with respect to time, 
and then multiplying the result by the square of the perpen- 
dicular wavenumber, k± = ak/R(z)- The average frequency 
Qk of the waves can be defined as an average over the power 



spectrum: 



LbPk{Co,z)dCb 
Pk(Co,z)dCb 



(57) 



We further average these frequencies over modes k to obtain 
the average frequency (1j„{z) for each bin n in wavenumber 
space. The three curves in Figure |8]; show uj„ for three dif- 
ferent heights in the reference model. Note that these average 
frequencies generally increase with perpendicular wavenum- 
ber, as expected for turbulent flows. The diamond symbol 
in Figure [8}; indicates the average frequency of the footpoint 
motions. 

In fully developed Alfvenic turbulence, the fluctuation are 
expected to reach a "critical balance" in which the aver- 
age wave period is of the order of the nonlinear transfer 
time, (^_lv_l)~', where vj, i s the average velocity. Following 
iGoldreich & Sridhari (IT995I) . we define a nonlinearity param- 
eter 

C„^^^i^, (58) 
where „ is the average velocity based on the power in bin 



n, which is given by 



Pv,„a±^„/Aa. Figure |8}l shows 



(„ as function of dimensionless perpendicular wavenumber. 
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Note that for low wavenumbers ^ 1, indicating the pres- 
ence of strong turbulence at all heights. At large wavenum- 
bers ^„ drops below unity, which is due to wave damping. 

These figures demonstrate that the waves injected into the 
corona through the TR have a broad range of perpendicular 
wavenumbers, and have temporal frequencies that are very 
different from the driver modes launched at the base of the 
photosphere. Therefore, the dynamics of waves and turbu- 
lence in the lower atmosphere is very important for under- 
standing the coronal heating problem. 

4.3. Open Field Model 

To demonstrate that the turbulence in the lower atmosphere 
is nearly independent of processes in the corona, we construct 
an open field model with the tube connected to the photo- 
sphere only at one end. In this model there is only a single 
TR, and the coronal part of the tube has constant temperature, 
density and field strength (Tq =1.5 MK, pQ = 9x 10"'^ g cm"^. 
Bear = 50 G). At the coronal end of the tube we use open 
boundary conditions such that the Alfven waves propagate out 
of the computational domain and no waves are injected. In the 
coronal part of this model there are only outward propagating 
waves, hence no nonlinear wave-wave interactions. However, 
wave reflection still occurs in the chromosphere and TR, so 
there are counter-propagating waves in the lower atmosphere, 
which generates strong turbulence. In Figure |9] various quan- 
tities are plotted as function of position along the flux tube 
for this open field model. Comparison of Figures |5] and |9] 
indicates that the energy density, heating rate and other pa- 
rameters of the turbulence in the lower atmosphere (r < 81.3 
s) are very similar to those in the closed field model. Figure 
|9}l shows the corona (region with r > 81.3 s) is still being 
heated, but this residual heating is due to the damping of the 
outward traveling waves, and is not due to turbulence. 

In the above simplified model for the open field, there is 
no turbulent decay of the Alfven waves in the corona. How- 
ever, more realistic models of coronal holes and other open- 
field structures have shown that there is significant reflec- 
tion of waves in the corona, and that Alfven wave turbu- 
lence can explain the heating and acceleration of the solar 
wind ( e.g., JD mitruk et al. 2002; Cranmer & van Ballegooijep 



wmg ( e.g.,i.lJ 
2001 ICrani 



ICranmer et al.l 120071: IChandran & HoUwegl l2009l 



Verdini et al.ll20Toh . 



4.4. Dependence of Heating Rate on Model Parameters 

We construct a series of closed field models with different 
values of the model parameters. Table [T] shows the subset of 
models in which the properties of the footpoint motions or the 
coronal loop length are varied. The first five columns in Ta- 
ble [T] show the model name, the correlation time tq, vorticity 
ojo and velocity Av',„,5 of the footpoint motions, and the loop 
length Lcor- Model 21 is the reference model described in the 
previous subsection, and models 22 - 27 are variants in which 
one of the parameters tq, wo or Lcr is changed relative to the 
reference model. For all models listed in Table [T] the coronal 
field strength Bcor = 50 G, the height of the transition region 
Ztr =1.8 Mm, the coronal pressure pc„r = 1.833 dyne cm"^, 
the time step Afo = 0.746 s, the expansion factor F = 1, and 
damping rate v^ax = 0.7 s"' . 

For each model we simulate the dynamics of the Alfven 
waves generated by random footpoint motions for a period of 
3000 seconds, and we compute the heating rate Q(z) averaged 
over the time interval t = [800,3000] seconds of the simula- 
tion. The period before t = 800 [s] is omitted because it some- 



times contains a large spike in heating that may be an artifact 
of the initial conditions. The heating rate Q(z) decreases with 
height in the lower atmosphere. The average heating rate in 
the chromosphere is determined as follows. The function Q(z) 
in the height range z = [700, 1300] km is fit with the following 
expression: 

Pq(z) 



chrom 



Pchro] 



(59) 



where pchmm is the density at z= \QQQkm.(pchmm ~ 3.6 x 10~" 
g cm"^), and Qcimm and rj are constants, which are determined 
by the fit. Therefore, Qchmm is the average heating rate at 
height z = 1000 km. We also measure the average coronal 
heating rate, Qcor- The values of 77, Qchmm and Qco,- are listed 
in the last three columns of Table [1] Based on the results in 
Table [T] the coronal heating rate can be approximated as 



Qa 



'2.91 X 10"' 0.45 + 
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[erg cm s ] , 



(60) 



where to is in seconds. Therefore, the heating rate increases 
nonlinearly with the vorticity luq of the footpoint motions, and 
decreases approximately inversely with loop length Lc„r. 

Table |2] shows the subset of models in which the coronal 
field strength and plasma pressure are varied. The pressure 
is controlled by the TR height, ztr (see section [J!2] i. The first 
six columns in Table[2]show the model name, the coronal field 
strength Bcr, the TR height ztr, the coronal pressure pcor, the 
loop length Lcor, and the time step Afo used in the simula- 
tion. Lcor is listed here because it varies slightly between these 
models (see Appendix B). The values of coronal field strength 
typically found in active regions lie in the range 10 - 500 G. 
However, the present version of the RMHD code has difficulty 
simulating Alfven waves in loops with Bcor > 200 G, which is 
due to the large Alfven speeds involved (v^ > 7000 km s"'). 
Therefore, in the present paper we only consider loops with 
Bcor in the range 12 - 200 G. For all models shown in Table |2] 
the footpoint motions are characterized by tq = 60 s, wo = 0.04 
s"\ and Avrms = 1-48 km/s. As before, the expansion factor 
F = 1 and the damping rate iy„,ax = 0.7 s"' . For each model we 
compute the heating rate Q(z) averaged over the time interval 
t = [800,3000] seconds of the simulation, and we derive the 
parameters 77, Qchrom and Qcor (see last three columns of Table 

Figures [TOk and [TOb show the coronal and chromospheric 
heating rates as function of coronal pressure pcor for six val- 
ues of coronal field strength (Bcor = 12, 25, 50, 100, 150 and 
200 G). The symbols show the simulation data listed in Ta- 
bled and the dashed curves show quadratic fits to these data. 
Both Qcor and Qchrom increase with coronal field strength. The 
decrease of Qcor for low coronal pressures may be explained 
as effects of wave reflection: for small pcor or high Bcor the 
coronal Alfven speed is very high, resulting in strong wave 
reflection in the chromosphere and TR. As a result, a larger 
fraction of the wave energy is dissipated in the lower atmo- 
sphere. Figure [Tob shows that the chromospheric heating rate 
depends only weakly on coronal pressure. 

As discussed in section [J!2l the peak temperature T„,ax along 
a coronal loop was chosen to be consistent with the first RTV 
scaling law, equation ([T]), but the heating rates Qcor found in 
the simulations are not necessarily consistent with the second 
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TABLE 1 

Dependence of Heating Rates on Footpoint Motions and Coronal 
Loop Length 



Model 


TO 








V 


Qchrom 


Qcor 




[s] 


[s-'] 


[km/s] 


[Mm] 




[erg/s/cm^] 


[erg/s/cm'] 


21 


60 


0.04 


1.48 


49.6 


0.346 


6.08 X 10-1 


2.97 X 10-3 


22 


100 


0.04 


1.48 


49.6 


0.399 


6.30 X 10-' 


2.30 X 10-3 


23 


200 


0.04 


1.48 


49.6 


0.440 


5.03 X 10-' 


1.89 X 10-3 


24 


60 


0.02 


0.74 


49.6 


0.434 


1.29 X 10-' 


0.90 X 10-3 


25 


60 


0.06 


2.21 


49.6 


0.350 


1.59 X 10" 


5.46 X 10-3 


26 


60 


0.04 


1.48 


25.6 


0.362 


6.01 X 10-' 


5.23 X 10-3 


27 


60 


0.04 


1.48 


99.6 


0.360 


6.47 X 10-' 


1.49 X 10-3 



„ 10^ 




: (b) 
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Fig. 9. — Model for Alfven wave turbulence in an open field. Various quantities are plotted as function of position along the flux tube: (a) kinetic and magnetic 
energy densities, and their sum Ed); (b) rms velocity Av,-,,,,,; (c) rms vorticity uinns', (d) kinetic and magnetic heating rates, and their sum Q{z)\ (e) rms magnetic 
fluctuation A6,-„,j; (f) rms twist parameter arms- Positions are given in terms of the Alfven travel time t(z). The corona is located in the range 81.3<r<118.6 
s. The TR is indicated by the vertical dotted line. 



scaling law, equation (|2]). Therefore, the models listed in Ta- 
bles [1] and |2] are not necessarily in thermal equilibrium. The 
solid curve in Figure [TOb shows the heating rate Qcor derived 
from equation dU with Lcor = 50 Mm. This curve shows the 
value of Qcor at which the coronal loop is in thermal equilib- 
rium (i.e., the heating is balanced by radiative and conductive 
losses). To the left of this curve, the time-averaged heating 
rate due to Alfven wave turbulence is larger than the radia- 
tive and conductive losses, so the coronal temperature will 
increase, more mass will be evaporated into the corona from 
the chromosphere, and the coronal pressure will increase. The 
opposite happens to the right of the solid curve. Therefore, al- 
though the present models do not include the effects of chro- 
mospheric evaporation, it is clear that in models with time- 
dependent coronal pressure the corona will have a tendency to 
approach the equilibrium state represented by the RTV scal- 



ing laws. 

The field-strength dependence of the equilibrium heating 
rate Qcor can be determined by finding the intersection points 
of the solid and dashed curves in Figure [TOh . and the corre- 
sponding chromospheric heating rate Qchrom can be obtained 
from Figure [TOb. The "equilibrium" heating rates are plotted 
in Figures [Tot andfTOh as function of coronal field strength, 
Bcor- Note that both heating rates increase with coronal field 
strength, and show no sign of saturation for large Bcor- The 
dependences on Bcor can be roughly approximated as power 
laws: 



2™, w 2.88 X 10" 



Qchrom ~ 6. 49 xlQ- 



Bcor 

50G 

1 / Bcor 

V50G 



0.55 



0.47 



[ergcm-^-'], (61) 



[erg cm s '], (62) 
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TABLE 2 

Dependence of Heating Rates on Coronal Field Strength and Plasma Pressure 



Model 




ZTR 


Pcor 




Ato 




Qcbrom 


Qcor 




[G] 


[Mm] 


r J / 9 n 

[dyne/cm J 


[Mm] 


[s] 




[erg/s/cm'j 


[erg/s/cm ] 


54 


200 


1.5 


4.585 


50.3 


0.296 


0.544 


1.268 


6.86 X 10"^ 


52 


200 


1.6 


3.375 


50.4 


0.252 


0.514 


1.357 


6.76 X 10"^ 


53 


200 


1.7 


2.487 


50.6 


0.216 


0.496 


1.386 


6.47 X 10"^ 


51 


200 


1.8 


1.833 


49.4 


0.186 


0.544 


1.236 


5.43 X 10-3 


55 


150 


1.5 


4.585 


50.0 


0.392 


0.501 


1.062 


5.07 X 10-3 


56 


150 


1.6 


3.375 


50.5 


0.337 


0.475 


1.144 


5.12 X 10-3 


57 


150 


1.7 


2.487 


49.4 


0.289 


0.489 


1.099 


4.84 X 10"3 


58 


150 


1.8 


1.833 


50.6 


0.247 


0.483 


1.017 


4.61 X 10"3 


38 


100 


1.6 


3.375 


50.6 


0.507 


0.427 


7.89 X 10"' 


3.75 X 10"3 


44 


100 


1.7 


2.487 


50.7 


0.433 


0.425 


8.05 X 10"' 


3.76 X 10"3 


29 


100 


1.8 


1.833 


49.5 


0.373 


0.436 


8.19 X 10"' 


3.42 X 10"3 


42 


100 


1.9 


1.351 


49.5 


0.319 


0.440 


8.00 X 10"' 


3.36 X 10"3 


40 


100 


2.0 


0.996 


50.5 


0.274 


0.431 


8.06 X 10"' 


3.03 X 10"3 


36 


100 


2.1 


0.734 


50.1 


0.235 


0.442 


8.03 X 10"' 


2.66 X 10"3 


30 


50 


1.6 


3.375 


49.3 


1.016 


0.405 


5.63 X 10"' 


2.90 X 10"3 


32 


50 


1.7 


2.487 


50.5 


0.864 


0.392 


5.99 X 10"' 


2.93 X 10"3 


21 


50 


1.8 


1.833 


49.6 


0.746 


0.346 


6.08 X 10"' 


2.97 X 10"3 


33 


50 


1.9 


1.351 


50.5 


0.637 


0.392 


5.83 X 10"' 


2.71 X 10"3 


34 


50 


2.0 


0.996 


50.6 


0.549 


0.360 


6.03 X 10"' 


2.64 X 10"3 


31 


50 


2.1 


0.734 


50.1 


0.469 


0.365 


6.01 X 10"' 


2.35 X 10"3 


35 


50 


2.3 


0.399 


50.4 


0.347 


0.371 


5.75 X 10"' 


2.06 X 10"3 


39 


25 


1.6 


3.375 


50.4 


2.018 


0.423 


4.63 X 10"' 


2.23 X 10"3 


45 


25 


1.7 


2.487 


50.3 


1.721 


0.403 


4.44 X 10"' 


2.13 X 10"3 


28 


25 


1.8 


1.833 


49.5 


1.481 


0.372 


4.61 X 10"' 


2.05 X 10-3 


43 


25 


1.9 


1.351 


50.3 


1.269 


0.390 


4.46 X 10"' 


2.05 X 10-3 


41 


25 


2.0 


0.996 


49.5 


1.100 


0.407 


4.45 X 10"' 


1.88 X 10-3 


37 


25 


2.1 


0.734 


50.4 


0.944 


0.422 


4.19 X 10"' 


1.85 X 10-3 


46 


25 


2.3 


0.399 


50.4 


0.694 


0.389 


4.31 X 10"' 


1.61 X 10"3 


50 


12 


1.9 


1.351 


49.6 


2.669 


0.414 


3.51 X 10"' 


1.39 X 10"3 


49 


12 


2.0 


0.996 


50.5 


2.283 


0.406 


3.70 X 10"' 


1.33 X 10"3 


48 


12 


2.1 


0.734 


50.4 


1.967 


0.382 


3.61 X 10"' 


1.30 X 10"3 


47 


12 


2.3 


0.399 


50.1 


1.438 


0.373 


3.52 X 10-1 


1.21 X 10-3 



although the fits are not very accurate at high field strength 
(see solid curves in Figures [Tob and [Toll). Combining equa- 
tions ( l60l l and dMT l. we obtain 



2.9 X 10" 



50 Mm) 



-0.92 



33 

0.45 -F — 

To 



1.65 



50 G 



1.48 km/s J 



[erg cm ^ s 



(63) 



where we replaced the vorticity by the velocity Av^s of 
the footpoint motions. Note that the equilibrium heating rate 
depends on coronal field strength Bcor, on loop length Lcor, 
and on the parameters of the footpoint motions (tq and Avrms)- 
The above expression assumes that the coronal temperature 
and pressure have the values predicted by the RTV scaling 
laws, and therefore neglects the e ffects o f the spatial and tem- 
poral variability of the heating (iKhmchuk. .2006.) . Also, it 
should be kept in mind that the numerical simulations have 
been done only in a limited range of coronal field strength 
and loop length (12 < B^or < 200 G, 25 < L„,r < 100 Mm), 
so the above expression for Qc,,,- is valid only within this lim- 
ited range. Furthermore, the above expression applies only to 
coronal loops with constant cross-section (F = 1), and the ef- 
fects of gravity have been neglected in the corona (but not in 
the lower atmosphere). Expression (|63) should not be applied 
to coronal loops on the quiet Sun, where magnetic fields are 
weaker and the effects of gravity and coronal loop expansion 
cannot be neglected. 



We also considered how the heating rates depend on the 
maximum damping rate Vma^ of the Alfven modes [see equa- 
tion (IB13I )1. For model Nos. 21, 28 and 29 we did a second 
simulation with twice the damping rate (I'max =1.4 s"'), and 
we found that the coronal heating rates Qc,,,- are changed by 
about -1-1%, -6% and -6%, respectively. The chromospheric 
heating rates Qchmm are changed by -6%, -1-2% and -10% for 
these three models. These numbers are relatively small com- 
pared to the factor 2 change in the damping rate, indicating 
that the rate of turbulent dissipation is insensitive to the value 
of the damping rate of the high wavenumber modes, as one 
would expect for a turbulent process. 

4.5. Effects of Coronal Loop Expansion 

Models of the coronal magnetic field (e.g., ISchriiver et al.l 
2005) predict that the field strength Bo generally decreases 
with height in the corona. Therefore, we expect that the 
cross-sectional area A of a coronal loop increases with 
height: A = $/Bo, where $ is the magnetic flux, which 
presumably is constant along the loop. However, observed 
X-ray and EU V loo ps o ften sh ow a pproximately constant 
cross-section jKlimchuk' llOOOl: IWaflto & KUmchukI [2000l : 
iLopez Fuentes et al. 2006i) . The reasons for this discrepancy 
are not well understood. The present paper cannot directly 
address this issue because our model does not include inter- 
actions between neighboring flux tubes, which we believe is 
important for resolving this issue. However, we can investi- 
gate how Alfven waves are affected by the expansion of the 
loop with height. We define the areal expansion factor F as 
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Fig. 10. — Dependence of coronal and chromospheric heating rates on coronal pressure (pcor) and field strength {Bcr) for a loop with constant cross-section 
and length Lcor = 50 Mm. (a) Average coronal heating rate Qcor as function of pcor- The symbols show results from the numerical models listed in Table|2]for 
Bcoi = 12 G iplusses), 25 G (crosses), 50 (stars), 100 G (diamonds), 150 G (triangles), 200 G (squares). The dashed curves show quadratic fits to these data, 
and the solid curve shows the heating rate predicted by the RTV scaling law, equation j2). The latter represents the condition of thermal equilibrium (balance 
between heating and cooling processes), (b) Heating rate Qchmm ^t height 7, = I Mm in the chromosphere, plotted as function of coronal pressure. Panels (c) and 
(d) show the equilibrium heating rates Qcor and Qduom as function of coronal field strength. The plus sign s sh ow th e va lues derived from the intersections of the 
dashed and solid curves in panel (a), and the solid curves show power-law fits to these data, see equations j6U and (62). 



TABLE 3 

Effect of Coronal Loop Expansion 



r 




V 


Qchyolu 


Qcor 


Fa(ztr) 




Pcor 


fcor 




[cm-3] 




[erg/s/cm^''] 


[erg/s/cm^] 


[erg/s/cm-] 


[erg/s] 


[erg/s] 




1 


0.44 X 1026 


0.346 


6.08 X 10-1 


2.97 X 10-3 


0.77 X 10' 


1.69 X lO^-* 


1.30 X 10" 


0.077 


2 


0.73 X 10^6 


0.386 


5.70 X ur' 


2.54 X 10-3 


1.09 X 10' 


1.66 X lO-"* 


1.85 X 10-3 


0.111 


3 


1.03 X lO^** 


0.353 


5.57 X ur' 


2.11 X 10-3 


1.25 X 10' 


1.66 X lO^"* 


2.17 X 10^3 


0.130 


4 


1.32 X lO^** 


0.377 


5.75 X 10-' 


1.76 X 10-3 


1.34 X 10' 


1.70 X 10^" 


2.32 X 10^3 


0.137 


6 


1.91 X 10^6 


0.371 


5.53 X 10-' 


1.43 X 10-3 


1.56 X 10' 


1.69 X 10^4 


2.73 X 10^3 


0.161 
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described in equation (ISTT i. This factor affects only the coro- 
nal part of the loop, and has no effect on the magnetic field 
Bo(z) in the lower atmosphere. We construct a series of mod- 
els with different expansion factors, but otherwise identical 
model parameters {ztr = 1.8 Mm, Lcr = 50 Mm, Bcr = 50 
G). For each model we compute the heating rate per unit vol- 
ume, Q(z), averaged over the time interval from 800 to 3000 
seconds and over the cross-sectional area of the flux tube. We 
also compute the rate of energy input into the corona. 



P^or=MZTRl)FA(ZTRl)-A(ZTR2)FA(ZTR2), 



(64) 



where A(z) = nR^iz) is the cross-sectional area of the tube, 
Ztr\ and ztr2 are the positions of the TRs, and Fa{z) is the 
time-averaged Alfven wave energy flux [see equation (IC12H . 
Finally, we compute two energy dissipation rates, P,oi and 
Pcor, which represent the dissipated power integrated over the 
total volume of the tube and over the coronal volume, respec- 
tively: 

rL rZTRi 

P,o, = / Qiz)A(z)dz, Pcor = / Q{z)A{z)dz. (65) 



-'0 'J Z.TR\ 

We find that Pcor and P'^^^ are equal to within a few percent, as 
expected for a statistically stationary state. This equality also 
shows that energy is very well conserved in these numerical 
models. 

The results of the simulations are summarized in Table |3] 
The nine columns of this table show the areal expansion fac- 
tor r, the coronal volume Vcor, the parameters rj and Qchmm 
describing the heating in the lower atmosphere [see equa- 
tion (|59]l1, the volume-averaged coronal heating rate Qcor, the 
Alfven wave energy flux Fa{ztr) at the TR, the integrated dis- 
sipation rates P,„, and Pcor, and the fraction of energy dissi- 
pated in the corona, fcor = Pcor/Ptot- The model with F = 1 
is the reference model discussed in section 14.11 We find that 
for the reference model only 7.7% of the available energy is 
dissipated in the corona; the remainder is dissipated in the 
photosphere and chromosphere at the two ends of the coronal 
loop. The energy flux Fa{ztr) increases with the expansion 
factor, and so does the coronal power Pcor and the fraction 
fcor- Specifically, for F = 6 about 16% of the available energy 
is dissipated in the corona. In contrast, the total dissipation 
rate P,o, is almost independent of the expansion factor. There- 
fore, the heating of the lower atmosphere is apparently not 
sensitive to the properties of the coronal part of the loop. 

The results of Table [3] show that, as F increases, more en- 
ergy enters into the corona and is dissipated at coronal heights. 
We suggest this increase is due to the fact that for large F the 
Alfven speed at the midpoint of the coronal loop (z = L/2) is 
drastically reduced compared to the case F = 1 . As a result, it 
takes a longer time for the Alfven waves to travel through the 
corona, and this gives the waves more time to dissipate their 
energy in the corona via turbulent processes. The wave energy 
is distributed over a coronal volume Vcor that increases more 
rapidly with F than the energy dissipation rate. Therefore, the 
average heating rate Qcor in the corona decreases with increas- 
ing F (see Table|3]l. 

Figure[TT]shows the magnetic field strength Bo(z) and time- 
averaged volumetric heating rate Q{z) as functions of position 
along the coronal loop. Note that Q{z) is relatively constant 
within the corona, and varies more gradually than the mag- 
netic field strength Bo(z). A power-law fit to these data for 
F = 4 and 6 yields Q{z) oc [Bo(z)]*' '*, so the exponent is sig- 
nificantly less than unity. We suggest this relative constancy 



of the volumetric heating rate may be due to Alfven wave 
reflection within the corona. For large F, the Alfven speed 
va(z) has a local minimum at the mid-point of the loop (z = 27 
Mm), and increases toward the ends. Hence, as the Alfven 
waves travel from the midpoint to the ends of the loop, the 
waves are partially reflected before they reach the TR. This 
enhances the wave energy density in the central part of the 
loop, and gives the wave more time to be dissipated near the 
midpoint. Clearly, wave reflections in the chromosphere, at 
the TR, and within the corona play an important role in de- 
termining the spatial distribution of the average heating rate 
Q(z). 

It is worthwhile to compare the computed heating rates Q(z) 
wit h those predicted from phen omenolo gical turbulence mod- 
els (Hossain et al. '1995'; Matthaeus et all 1999 ": Dmitru k et"!!] 
2001, 2002; Chandran et al. 2009). In such models the fluc- 
tuations are expressed in terms of the Elsasser variables, here 
defined as Z± = \ zLB±/ \/4tt po. The turbulent dissipation 
rate is given by equation (57) of lCranmer & van BallegooijenI 
(l2005h : 



Qpheniz) = PO 



z^z++zlz+ 



(66) 



where L±(z) is the outer scale of the turbulence, and the quan- 
tities Z±(z) are the rms values of the Elsasser variables. The 
latter are averaged over the cross-section of the flux tube [see 
equation ( ICBI ll and over time. We compute Z±{z) for sev- 
eral models and find that the turbulence is significantly un- 
balanced (Z+ y Z_), which is due to the strong dissipation in 
our models. The perpendicular length scale is set equal to 
the flux tube radius, L± = R{z)- Figure [12] shows the ratio 
q{z) = Q{z) / Qpheniz) plotted as function of position along the 
flux tube for three values of the coronal loop expansion factor 
(F = 1, 3 and 6). Positions are given in terms of the Alfven 
travel time t{z) from the left footpoint. We find that q w 0.2 
in the photosphere and low chromosphere, and the ratio peaks 
in the upper chromosphere, indicating the turbulence is more 
efficient there. For the reference model q w 0.3 in the corona 

i solid cuTwe), similar to the value of 0.5 used by iBreech et al.l 
20091) in their model of the solar wind. Similar values for 
q are found for the models with larger expansion factors (the 
downward spikes just above the TR may be due to insuffi- 
cient spatial resolution). We conclude that q is significantly 
less than unity at all heights. 

5. SUMMARY AND DISCUSSION 

According to the present model, the interactions of turbu- 
lent convective flows with kilogauss flux tubes in the pho- 
tosphere produce transverse displacements of magnetic field 
lines on length scales less than the width of the flux tubes, 
i.e., less than about 100 km (see FigurelTJ. We assumed that 
the photospheric footpoint motions have velocity amplitudes 
in the range 1-2 km/s and correlation times of 60 - 200 s. 
We studied the dynamics of the Alfven waves produced by 
such footpoint motions, and found that such motions produce 
Alfvenic turbulence at larger heights in the flux tube. The 
predicted dissipation rates of the turbulence as function of 
height are sufficient to reproduce the observed rates of chro- 
mospheric and coronal heating in active regions. We conclude 
that fine scale magnetic braiding can drive Alfvenic turbu- 
lence and can produce sufficient heating for both the chro- 
mosphere and corona. However, we should emphasize that 
we do not have any direct observational evidence for the exis- 
tence of the assumed footpoint motions, as scales of 100 km 
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Fig. 1 1. — Effects of coronal loop expansion, (a) Magnetic field strength Bq(z), and (b) volumetric heating rate Q(z), as functions of position z along the 
coronal loop for different values of the loop expansion factor: F = 1 (full), F = 2 (dotted), F = 3 (dashed), F = 4 (dash-dotted), and F = 6 (dash-triple-dotted). 
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Fig. 12. — Ratio of the numerically computed heating rate Q(z) and the 
rate Qphen(i) predicted by phenomenological turbulence models. The ratio 
q = QlQphen is plotted as function of Alfven travel time t(z) from the left 
footpoint for three different values of the loop expansion factor: F = 1 (full), 
F = 3 (dotted), F = 6 (dashed). Only half of the loop is shown (0 < z < L/2). 

or less cannot be resolved with present-day solar telescopes. 
Therefore, our conclusion can be made only for the chosen 
values of the photospheric driver; for higher or lower values 
the heating does not meet the observational constraints. Our 
results indirectly provide a constraint on the amplitude of the 
footpoint motions necessary for the proposed mechanism to 
heat the chromosphere and corona. Our model produces a 
corona l heating profile that is similar to that of a nanoflare 
storm (lKlimchukll2006h in the spatial and temporal distribu- 
tion of heating events, and in the final source of the energy 
released - dissipation of magnetic fields. The difference is 
that we have provided a fully self-consistent source for the 
nanoflare energy. The source depends critically on the exis- 
tence of a chromospheric reservoir of turbulent Alfven waves. 



To investigate the nonlinear dynamics of the Alfven waves, 
we considered a flux tube that extends from the photosphere 
through the chromosphere into the corona. We mainly con- 
sidered closed coronal loops, and for simplicity we assumed 
that a single kilogauss flux tube at one end of the loop is con- 
nected to a single flux tube at the other end. Therefore, we 
ignored the fact that on the real Sun the photospheric flux 
concentrations at the two ends are uncorrected and do not 
perfectly match up. Furthermore, we assumed that the flux 
elements do not split up or merge during the simulation, so 
that the flux tube retains its identity. The pressure, density 
and magnetic field strength were taken to be fixed functions 
of position z along the flux tube. The assum ed structure of 
the lower atmosphere is based on model P by Fontenl a et al.l 
(|2()09); this model represents faculae in active regions (i.e., 
very bright plage). The chromosphere-corona TR was treated 
as a discontinuity where waves can reflect, and we neglected 
the gravitational stratification in the corona, so the pressure is 
constant along the coronal part of the loop. 

The dynamics of the waves and their structure in the 
plane perpendicular to t he loop axis we re described using the 
RMHD approximation dStraussI [l976b . which assumes that 
the magnetic and velocity perturbations of the waves are per- 
pendicular to the mean magnetic field of the flux tube. We 
found that the Alfven waves strongly reflect as they propa- 
gate through the chromosphere. This reflection is due to the 
increase of Alfven speed with height. Even stronger reflec- 
tion occurs at the TR where the Alfven speed suddenly in- 
creases by about a factor of 15. These reflections produce in 
the chromosphere a pattern of counter-propagating waves that 
are subject to nonlinear wave-wave interactions. We found 
that the waves quickly decay into turbulence, causing most 
of the wave energy to be deposited in the lower atmosphere 
(photosphere and chromosphere). Only a small fraction of 
the wave energy is transmitted into the corona. Energy is 
injected at both ends of the coronal loop, producing in the 
corona two sets of counter-propagating waves that reflect off 
the TRs and decay into turbulence. The dissipation is en- 
tirely due to Alfven wave turbulence, and is not due to phase 
mixin g or re sonant absorption (e.g., Hevvaerts & P i'ie stj |1983l; 
iDe Groof & GoossensH2002l) . which are not included in the 
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present model. 

We constructed a variety of loop models with different val- 
ues of the model parameters. For the reference model (section 
14.1b we assumed the photospheric footpoint motions have an 
rms velocity of 1.48 km/s and a correlation time tq = 60 s. We 
found that the heating rate Q(z, t) varies strongly with position 
z along the loop and with time, which is due to bursts of wave 
activity and subsequent dissipation of the waves via turbu- 
lence. Only a small fraction of the wave energy is dissipated 
in the corona (about 7.7%), and the remainder is dissipated in 
the lower atmospheres at the two ends of the modeled loop. 
We computed the time-averaged heating rate Q{z) and found 
that this rate decreases with height in the lower atmosphere 
from about 10 erg cm"-' s"' in the photosphere to about 0.2 
erg cm~^ s~' in the upper chromosphere. These values can 
be compared with rates of radiative loss found in semiem- 
pirical models of the so l ar atmosphere (e.g., Avrett 1985J 
Anderson & Athav'll 989t iFontenla et al.lfl999l llooa ,20091) . 



Fontenlaet al. (2009) computed radiative loss rates as func- 
tion of gas pressure (i.e., height) for several different wave- 
length bands (see their Figure 14). For model P the largest ra- 
diative loss occurs in the 2000-3000 A band, which contains 
the Mg II resonance lines. The total losses over all wave- 
lengths are in the range 1-2 erg cm"^ s"\ somewhat larger 
than the values found in the present model. Therefore, our 
model may underestimate the amplitude of the Alfven waves 
and rate of energy dissipation in the chromosphere. 

In the corona the reference model predicts Qcr ^ 3 x lO""* 
erg cm"^ s~', which is sufficient to explain the heating of typ- 
ical active region loops. The predicted velocity amplitude of 
the Alfven waves in the corona is in the range 20 - 40 km/s, 
similar to the values fou nd in observations of spectral Une 
widths in active region s (iDere & Masonl[T993l: IWarren et al.l 
l2008tlLi & Ding||2009l) . We investigated how the coronal and 
chromospheric heating rates depend on the model parameters 
(see section |4~4] |. We used the RTV scaling laws to determine 
the "equilibrium" heating rates for which the coronal heating 
is balanced by radiative and conductive losses, and derived a 
scaling law describing how the coronal heating rate Qcor de- 
pends on the coronal field strength, loop length, and the pa- 
rameters of the footpoint motions [equation (l6?tl. This result 
may be compared w i th other models of co ronal heating (see 
iMandrini et ani2000t[Schrijver et al.ll2004 . 

We also considered coronal loops with non-constant cross- 
section (section 14.5b and found that both the energy flux 
Fa{ztr) entering the corona and the total power Pc„r dissipated 
in the corona increases with the areal expansion factor F. The 
increase of Per with F is due to the decrease in the coronal 
Alfven speed, which lengthens the Alfven travel time in the 
corona and gives the waves more time to dissipate their en- 
ergy via turbulence. The fraction of energy dissipated in the 
corona also increases with F, but is still only about 16% for 
the model with F = 6. The average coronal heating rate Qcor 
decreases with F, which is due to the fact that the coronal vol- 
ume increases faster than the dissipated power Pcor- We found 
that for models with F > 1 the heating rate Q{z) varies with 
position along the coronal part of the loop, and decreases to- 
wards the midpoint of the loop where the field strength Boiz) 
is lowest. However, Q{z) varies less than Bo{z), so the heating 
is not strongly concentrated near the coronal base. This ef- 
fect can be understood in terms of reflection of the waves by 
gradients of Alfven speed in the corona. 

Observations of the active corona using the Transition Re- 



gion and Coronal Explorer (TRACE, see iHandv et al.lfT999h 
have shown that the corona is highly dynamic and full of 
flows and wave phenomena (e.g..'S chriiver et aLlll999l) . This 
has led some authors to suggest that the heating of active re- 
gion loops is localized in the low corona and involves ener- 
gization proce sses that operate in t he chromosphere and tran- 
sition region ([Aschwandenl 120011: [Aschwanden et al. 200'!^ 
iDe Pontieu et all l2009t) . However, .KUmchuk et al., t20ldh 
show that heating concentrated in the low corona would lead 
to thermal non-equilibrium of the plasma, and such non- 
equilibrium is not consistent with the observed properties of 
warm (1-2 MK) loops observed in active regions. Instead, 
they suggest coronal loops are heated impulsively by storms 
of nano flares that are distributed throughout the coronal loop 
(also see |Parker"1988| |Lu_&Hamilton 19911 lKlimchukll2006l 
l2009l ). The present model has some features in common with 
both types of models. On the one hand, we find that the aver- 
age heating in the chromosphere is 2 or 3 orders of magnitude 
larger than that in the low corona, and that most of the heating 
occurs in the lower atmosphere, which is consistent with the 
ideas of De Pontieu et al. (2009). On the other hand, in the 
coronal part of the loop the heating depends only weakly on 
position, so the lo op is likely to be thermally stable, consistent 
with the results of lKlimchuk et alj (^010). 

Although the present model was described in terms 
of the propagation and dissipation of Alfven waves, the 
waves reflect at various heights and undergo nonlinear 
interactions, producing twisted and braided fields similar to 
those f ound in previous bra i ding m o dels (e.g., [M ikic et aj] 
'1989"; 'Lonscooe & Strauss 1994; Hendrixetab 1 1996 : 
Galsgaard & Nordlund 1996: Ng & Bhattacharieel Il998 : 
ICraig & Snevdl l2005t iRappazzo et al.1 120081) . Also, flie 
model includes the effects of magnetic reconnection, and 
the time-dependent heating seen in Figure |4}l can be thought 
of as a sequence of "nanoflares" (Parken il988l) . However, 
unlike in previous "quasi-static" braiding models the twisted 
structures are highly dynamic and are far from the force-free 
equilibrium state. For example, the twist parameter a is not 
at all constant along the field lines (see Figure |5f), as would 
be expected for a force-free state. This is due to the inclusion 
of the lower atmospheres at the two ends of the coronal loop, 
which produces turbulence in the chromosphere. Also, the 
transverse motions of field lines at the photospheric base 
occur entirely inside the magnetic flux elements, hence these 
motions do not contribute to the "random walk" of the flux 
elements on the photosphere. This has the advantage that the 
velocity amplitudes can be relatively large, yet be consistent 
with obs ervational constrain t s on the photosphe ric diffusion 
constant (iDeVore et al.llT98llBerger et al.l[T998h . 

According to the present model, the transverse scale of the 
braids in the corona is very small, £x < Rcor = 529 km, so 
the braids would likely not be observable with existing X- 
ray or EUV telescopes. Furthermore, the rms value of the 
transverse field fluctuations is small compared to the mean 
field, ABrms/Bo ~ 0.025, so the mis-alignment angle between 
the field lines and the tube axis is only about 1.4°. Both the 
transverse scale of the braids and the mis-alignment angle are 
consistent with the fact that no braiding is observed in the 
corona on scales of several arcseconds or larger (see section 
12.1b . The mis-aUgnment angle is much smaller than that re- 
quired in qua si-static braid i ng models th at match the observed 
heating rate ('Parker'19831 iPriest et al.l 12002). We conclude 
that magnetic energy is injected into the corona via a dynamic 
(not quasi-static) braiding process. The existence of turbu- 
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lence in coronal loops may have importan t implications for 
the d amping of coronal loop osci l lations dNakariakov et alJ 
ll999HOfman & Aschwandenl2 002':'Roberts'2008') and for the 
cross-field diffusion of electrons (Galloway et al. 2006). 

In the present paper we did not consider the plasma re- 
sponse to the heating. In future work we will construct mod- 
els of the thermal structure of coronal loops heated by Alfven 
wave turbulence, including the effects of thermal conduction 
and radiative losses. The results will be compared with pre- 
dictions from quasi-static braiding models. Another line of 
future research will be to investigate how the temporal fluctu- 
ations of the heating rate affect the temperature and density. 
Is thermal conduction strong enough to smooth out the spa- 
tial and temporal fluctuations? Ideally, such studies should 
take into account the non-uniformity of the temperature and 
density in the perpendicular plane; this may require a more 



complete 3D MHD model. We also intend to consider the 
interactions between neighboring flux tubes in the lower at- 
mosphere, including the effects of the splitting and merging 
of kilogauss flux tubes. 
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APPENDIX 

A. COMPARISON WITH TRANSPORT EQUATIONS BASED ON ELSASSER VARIABLES 

In section [TT] we formulated the Alfv en wave transp ort equations in terms of scalar fields /, h, uj and a. In contrast, previous 
works often used the Elsasser variables (lElsasserlll950b . which are vector fields: 



Z±(r,f) = Vi± 



B 



(Al) 



where \± = V/ x Bq is the velocity, and B^ = V/i x Bq is the perturbation of the magnetic field. To demonst rate that the present 
formu lation of the transport equations is equivalent to that used by others, we start with equation (21) of IZhou & MatthaeusI 
([T990h : 



= -(UtV^)-VZ±- 
-Z^- ( vu±- 



1 



-z^)v.(iu±v^) 



VBo VP-Z^-VZ±, 



PO 



(A2) 



where = Rp/^^/Anpg is the Alfven velocity vec t or, U is the outflow v elocity, and P is the p erturbation in the total pressure 
(also see lDmitruk et al.l[200lHVerdini & Wiill2007l: iBuchUn & V elli 2007: ICranmer et al.ll2007l) . In the present paper we neglect 



flows along the field lines, so we set U = 0. The terms in equation 
^ and (fT2l) . which yields 

dZ± ^ Bod 
— — = ±vaBo • VZ± T -IT 3- 
ot 2 dz 



involving Va and Bq can be simplified by using equations 



1 



1 



dBo 1 
2 V47r/9o dz po 



Z 



VZj 



(A3) 



We now take the curl of the above equation and evaluate the component parallel to the mean magnetic field. Using equation (fT2l l 
to evaluate the derivatives of Bq, we find 

duj± 



^^■ = ±v,(^Bo.Vc.±-|| 



BQd_ / 1 

2 dz \ y/4-npo 



± 



1 1 dBo 



-Wzr 



-Zzr-Vw±±7V, 



2 y/Anpo dz 

where w± = Bq • (V x Z±) are the vorticities associated with the Elsasser variables, and J\f is given by 



dZ- V dZ 



'+■)■ 



dZ-^y dZ+,y ^ dZ- 



9Z+..V ^ dZ-y dZ+_x 



dx dx dx dy dy dx dy dy 



(A4) 



(A5) 



Combining the linear terms in equation (IA4l i and using (lo^-lo-)/! = v/^a, we find that the above equations are identical to 
equations (|44] | and ( |45] l. Therefore, the present formalism is equivalent to that used in previous work based on Elsasser variables. 



B. NUMERICAL METHODS 

The four quantities w, a, h and / are determined by equations (l33T l. (l35t . (|40] | and (l39t derived in section [TTI In this Appendix 
we describe the numerical methods used for solving these equations. Let r be the radial distance to the tube axis (r < R{z)), ^ 
the azimuth angle, and z the distance along the tube. It is convenient to write the stream function as f{^,if,z,t), where ^ = is 
the fractional distance, and similar for h, uj and a. Then spatial derivatives along the background field can be written as partial 
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Fig. 13. — The 92 basis functions Fi(x,y) for Omax = 20. The two driver modes are shown as the seventh and eighth image of the top row. 



derivatives at constant for example, Bo • V/ = df /dz in equation (l39T l. In this way the vector Bo can be ehminated from the 
equations. 

The transverse structure of the waves is described using a spectral method. We use a finite set of basis functions Fi{^, ip), where 
index / enumerates the basis functions (/ = 1 , • • • ,A^). The functions are assumed to vanish at the tube wall, Fi(l, ip) = 0, consistent 
with our requirement that v and B are parallel to the wall. The basis functions are mutually orthogonal and are normalized such 

2 

7o Jo 

where 5ij is the Kronecker delta symbol. An arbitrary spatial distribution /(^, Lp) can always be projected onto the basis functions 
by evaluating 

N 
1=1 

where ft ai^e the mode ampUtudes: 

fi = - 1 I m, vWiiL (B3) 

Jo Jo 

In the following we assume / w /pioj, i.e., we suppress high spatial frequencies that are not present in the basis functions. The 
basis functions /^(^, (p) are chosen to be the eigenmodes of the operator; 

R^V\Fi = XiFi, (B4) 

where A,- is the dimensionless eigenvalue. The solutions f) are proportional to JnifliOi where m = nii is the azimuthal mode 
number (m > 0), Jm{x) is the Bessel function of the first kind, and a, is a zero of the Bessel function, Jmiad = 0. The modes are 
ordered according to their m-values (first the m = modes, then m = 1 modes, etc.), and for fixed m they are ordered according 
to their a-values. The non-axisymmetric modes always come in pairs, one proportional to cos(mp), the other proportional to 
sm{mp), which we count as two separate modes / and / + 1, respectively. The eigenvalues of the modes are A, = -aj. We find all 
modes with a, below a prescribed maximum value, a„M.v For example, for a,„a.x = 20 we find N = 92 modes (6 modes with m = Q, 
12 modes with m = 1, 10 modes with m = 2, etc., up to m = 15). The basis functions are shown in Figure [13] 
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As shown in section ITTl all nonlinear terms in the RMHD equations involve the bracket operator, defined in equation (HTt . In 
polar coordinates this quantity can be written as 

1 f df ds df dg\ 



r \ dr dip dr ^ 

where /(r, tp) and g(r, ip) are arbitrary functions. Expanding /, g and b in terms of basis functions, we obtain 

^ N N 

j=i i=i 

where Mkji is the coupling matrix: 

1 f^'' fdFjdFi dFidFi\ 

Note that Mi^jj is anti-symmetric with respect to / and j, and using partial integration one can show that M/tji = Mikj, where we use 
dFi/df = for ^ = 0. It follows that the six matrix elements that couple any three modes j and k are closely related: 

Mkji = Mikj = Mjik = -Mkij = -Mijk = -Mjki. (B8) 

To compute these elements we only need to consider the case k> j > i (the "unique" elements). The three modes must have 
compatible cos(m(p) or sm{mip) dependence, and for m, > we only need to consider the case = nij + mi. We find that for 
cimax = 20 there are 7662 unique matrix elements. 
The quantities /, h, to and a can be written in terms of basis functions, for example: 

N 

M,p,z,t) = Y,fi(z,t)Fi(^,v), (B9) 

i=I 

where fi{z,t) are the amplitudes of the different transverse modes. Inserting these expressions into equations (l40t and (|39t . and 
using expression (IB6I 1. we find 

^ n , N N 

du}k_ 2 00^k 

dt ""^ dz ' 



■ = -TT^ + ^ '^Mkjiivlajhi - ujjfi) - , (BIO) 



j=i i=i 

t = f44££""*'--*. <B.l) 

j=l 1=1 

where 

ak = iak/Rfh, and ujk = (ok/Rffk, (B12) 

and we added artificial damping terms. The damping rate i^k is assumed to be independent of time t and position z along the tube, 
but depends on the perpendicular wavenumber to the sixth power ("hyperdiffusion"): 

Vk = ^max\-^] , (B13) 

where Vmax is the damping rate of the highest wavenumber modes (we use Vmax = 0.7 s"'). The alternative form of the dynamical 
equations, given by equation (|44] |. yields 

duj±,k duj±,k dvA 

—— = ±va^ VA—ak+---. (B14) 

at oz dz 

Here we omit the nonlinear and damping terms because those terms will be evaluated directly from equations (IB 10) and (IB 1 11 1 
(see below). 

In order to accurately describe the Alfven wave propagation, we use a grid of positions z„ such that the Alfven travel time At 
between neighboring grid points is exactly equal to the time step Afo of the calculation: 

= Afo = constant. (B 1 5) 

VA.n+l/2 

Here ,,+1/2 is the average speed between grid points n and n + 1. The Alfven travel time from the TR to any position z in the 
corona is 

t(z)-t(ztr) = -^/[«(z)], (B16) 

VAiZTR) 

where u(z) = -l+2(z-ZTR)/Lcor, and f(u) is defined by 

r" / 1 — 8m^\ "'''^ 

M=J ( no ) [l + (T-m-u^)]du. (B17) 
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Here we used equations (ISOl l. dSTT i and ( l53T l, and va(z7-s) is the Alfven speed at the coronal base, va{_Ztr) = Bcor/ V^Trpofo/?)- 
The grid of u„ values is found by inverting the function f{u). To ensure that the total Alfven travel time through the corona is a 
multiple of Afo, we make a slight adjustment to the assumed coronal loop length Lcor- 

We now describe the procedure for solving equations (IBIO) and (IBllI) . Let ujk,„(t) be the vorticity for mode k at grid point 
n and time t, and similar for the magnetic flux function hk,„{t). To advance these quantities from time f to f' = f + Afo, we use 
operator splitting, i.e., each time step Afo consists of two parts. In Part 1 we compute the effects of nonlinear mode coupling 
and damping, and in Part 2 we determine the effects of wave propagation and reflection. For Part 1 we can omit the propagation 
and reflection terms, i.e., the first term on the right-hand side of equation ( IB 1 Oi l and the first and second terms in equation (IBl II) . 
The resulting equations no longer contain any z-derivatives, so the equations can be solved separately for each position n along 
the loop. For Part 1, advancing uj^j, and hk.,, from time t to t' is done using a fifth-order Runge-Kutta method (P ress et al.lll992i) . 
The method uses adaptive stepsize control, i.e., the time interval Afo may be broken up into two or more smaller intervals with 
Af < Afo. On the other hand, for Part 2 (wave propagation and reflection) we can omit the nonlinear terms in the equations, so 
we instead use the linearized form of equation (IB14l i. where the Elsasser-like variables a;± jt,„ are computed from equation (i43] l. 
The solution of equation (IB 141) for time step Afo can be written as 

'^±.t,n(f + Afo) = tJ±,(t.„±i(f)-Afo ( ) Q;i,„±i/2(f), (B18) 

"2 /«±l/2 



where a^ „±i/2 is the average value of a in between neighboring grid points. The first term in equation ( IB18I ) describes wave 
propagation, and the second term describes the coupling due to gradients in Alfven speed (gradients in the TRs are treated 
separately, see below). Note that step 2 can be done separately for each mode k because there is no mode coupling in the 
linearized equations. Therefore, in each time step Afo the effects of both linear and nonlinear terms in the equations are taken 
into account. 

In our model there are two chromosphere-corona TRs, one at each end of the coronal loop. At these TRs the mass density (and 
therefore the Alfven speed) changes discontinuously with position, which causes strong wave reflection. Let ztr be the position 
of one of these TRs. In our model there are two grid points associated with this position: z\ just below the discontinuity, and 
Z2 just above it. Let va.\ and va.i be the Alfven speeds at these points, and let fi ± and f2.± be the corresponding inward and 
outward propagating waves (note that i < v^.i for the first TR, and i > v^o for the second one). The waves /i + and /2,_ that 
propagate away from the TR are assumed to be linear combinations of the waves /; _ and /2 + that propagate towards it: 

/l,+ = C„/i,_ + Ci2/2,+ , (B19) 
/2,- = C2l/l,- + C22/2,+ , (B20) 

where the c,/s are reflection and transmission coefficients. Then the stream functions in the two regions are given by 

/i = 5[(1+cii)/k- + ci2/2J, (B21) 

/2 = i[C2l/l,- + (l+C22)/2,+], (B22) 

h = \[-{\-cxi)fi,- + cnfiAlvAA, (B23) 

hi = 5[-C2l/l.- + (1 -C22)/2,+ ]/vA.2. (B24) 

At Ztr both the velocity and the magnetic field perturbations must be continuous (see IChandrasekhad 1 96 ll) . so we require /i = fi 
and h\ = hi- Moreover, these conditions must be satisfied for any value of /i _ and /2.+. It follows that the coefficients c,y are 
given by 

c-1 2 2c 

Cll =-C22 = — — -, Ci2 = — — -, C2i = — — -, (B25) 
c+1 c+1 c+1 

where c = va.i/vaa is the ratio of Alfven speeds. 

The code was tested in various ways. We studied the propagation of Alfven wave packets in a uniform flux tube to verify that 
no wave dispersion occurs. We also considered nonuniform tubes and verified that energy is conserved in wave reflections. For 
the stratified loop we verified that the total wave energy (integrated over the entire volume of the tube) is constant in time when 
a wave packet is present in the initial conditions, the footpoints are held fixed, and the damping rate vq = 0. This demonstrates 
that wave energy is conserved by the nonlinear interactions. For the full model we verified that the energy injected by footpoint 
motion minus the energy dissipated by damping equals the rate of change of the total energy to an accuracy of about 5% of the 
net input rate. 

C. ENERGY EQUATION FOR ALFVEN WAVES 

We first derive expressions for the averages of vector quantities over the cross section of the flux tube. The rms velocity 
Avrm.5(z,f) is given by 

(Av™)2^< |v|2>=— / / \Wf\^ rdrd^=—^ / f lj rdrd^ 
T^R Jq Jo t^R Jo Jo 

. .2^ .1 / N \ / 1 ^ \ 1 ^ 

= - / / T.fjpj ^ E«'/'^' = ^ E«^v. 

JO JO y J \ / 



f, (CI) 
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where < • • • > denotes the spatial average over the cross section of the tube. Here we used partial integration to obtain w 
(= -V\f), and we used equations (IB II) . ( IB9I ) and (IB12I) . Similarly, the rms magnetic fluctuation ABrms(z,t) is given by 



^_ J 



and the rms values of the Elsasser variables are 



^<|z±p>= 



1 ^ 



k=l 



It foUows from equations ( ICll l and iC2\ that the contributions of mode k to the velocity and magnetic power spectra are 



PvAz, t) = (at/Rffi and Pb^z, t) = B'^{ak/Rfhi. 



- d2/ 



^2^.2 



(C2) 



(C3) 



(C4) 



We now derive the energy equation in the context of the present model. The kinetic and magnetic energy densities are given by 



Ekiniz, t) = 2 po( Av,„„) = ^ 2^ ajj^ 



N 



k=\ 



Multiplying equation (IB 1 Oi l by po/t and summing over k, we obtain the following equation for the kinetic energy density: 



dEkin 



k=\ 



Bl d 



,2"0 



N N 



kin ; 



where Qkin is the rate of kinetic energy loss due to damping: 



Qkm(z,t) 



El 

R^ 



(C5) 
(C6) 

(C7) 

(C8) 



k=l 



Similarly, multiplying equation (IBl lb by B\/(ATi)(ak/RYhk yields an equation for the magnetic energy density: 



dE,r. 



dt 



k=l 



'R2 Ait 



N N 



dfk I dBQ 1 V^V^ ,^ . , 



mag J 



where Qmag is the rate of magnetic energy loss: 



Adding equations (|C7| i and ( |C91 l yields 



Qnmg{z,t) = --^ ^ Vkulhl. 



AttR^ 



k=l 



dE d /Fa\ 



(C9) 



(CIO) 



(Cll) 



where E(z, t) = Eiii„+E,„a^ is the total energy density, Q{z, t) = Qiij„ + Qmag is the total dissipation rate, B(,{z) represent the diverging 
geometry of the background field, and FA{z,t) is the Alfven wave energy flux averaged over the cross-section of the flux tube: 



4ttR^ 



(C12) 



k=l 



Note that the terms involving Mtj, drop out of the energy equation (ICl It . This is a result of equation ( IB8b . which shows th at the 
matrix Mtj, is anti-symmetric with respect to interchange of any two indices. Therefore, the nonlinear terms in equations (IBIO) 
and (IBl lb do not have a direct effect on energy transport along the tube, but they are of course responsible for exciting the high 
wavenumber modes that are subject to damping (Q). Therefore, our numerical scheme using basis functions /^(^, ip) is excellent 
for studying the energetics of the waves. 
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